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Abstract 



We present a comprehensive semiclassical investigation of the three- 
dimensional Sinai billiard, addressing a few outstanding problems in "quan- 
tum chaos" . We were mainly concerned with the accuracy of the semiclas- 
sical trace formula in two and higher dimensions and its ability to explain 
the universal spectral statistics observed in quantized chaotic systems. For 
this purpose we developed an efficient KKR algorithm to compute an ex- 
tensive and accurate set of quantal eigenvalues. We also constructed a 
systematic method to compute millions of periodic orbits in a reasonable 
time. Introducing a proper measure for the semiclassical error and using 
the quantum and the classical databases for the Sinai billiards in two and 
three dimensions, we concluded that the semiclassical error (measured in 
units of the mean level spacing) is independent of the dimensionality, and 
diverges at most as log h. This is in contrast with previous estimates. The 
classical spectrum of lengths of periodic orbits was studied and shown to 
be correlated in a way which induces the expected (random matrix) corre- 
lations in the quantal spectrum, corroborating previous results obtained in 
systems in two dimensions. These and other subjects discussed in the re- 
port open the way to extending the semiclassical study to chaotic systems 
with more than two freedoms. 

PACS numbers: 05.45.+b, 03.65. Sq 

Keywords: Quantum chaos, billiards, semiclassical approximation, Gutzwiller 
trace formula. 
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1 Introduction 



The main goal of "quantum chaos" is to unravel the special features which charac- 
terize the quantum description of classically chaotic systems JT], 0. The simplest 
time independent systems which display classical chaos are two dimensional, and 
therefore most of the research in the field focused on systems in 2D. However, 
there are very good and fundamental reasons for extending the research to higher 
number of dimensions. The present paper reports on our study of a paradigmatic 
three-dimensional system: The 3D Sinai billiard. It is the first analysis of a sys- 
tem in 3D which was carried out in depth and detail comparable to the previous 
work on systems in 2D. 

The most compelling motivation for the study of systems in 3D is the lurking 
suspicion that the semiclassical trace formula || - the main tool for the theo- 
retical investigations of quantum chaos — fails for d > 2, where d is the number 
of freedoms. The grounds for this suspicion are the following || . The semiclassi- 
cal approximation for the propagator does not exactly satisfy the time-dependent 
Schrodinger equation, and the error is of order h 2 independently of the dimension- 
ality. The semiclassical energy spectrum, which is derived from the semiclassical 
propagator by a Fourier transform, is therefore expected to deviate by 0(h 2 ) from 
the exact spectrum. On the other hand, the mean spacing between adjacent en- 
ergy levels is proportional to h d || for systems in d dimensions. Hence, the figure 
of merit of the semiclassical approximation, which is the expected error expressed 
in units of the mean spacing, is 0(h 2 ~ d ), which diverges in the semiclassical limit 
h — > when d > 2! If this argument were true, it would have negated our abil- 
ity to generalize the large corpus of results obtained semiclassically, and checked 
for systems in 2D, to systems of higher dimensions. Amongst the primary vic- 
tims would be the semiclassical theory of spectral statistics, which attempts to 
explain the universal features of spectral statistics in chaotic systems and its rela- 
tion to random matrix theory (RMT) [Q, |5|. RMT predicts spectral correlations 
on the range of a single spacing, and it is not likely that a semiclassical theory 
which provides the spectrum with an uncertainty which exceeds this range, can 
be applicable or relevant. The available term by term generic corrections to the 
semiclassical trace formula @, [7|, |§ are not sufficient to provide a better estimate 
of the error in the semiclassically calculated energy spectrum. To assess the error, 
one should substitute the term by term corrections in the trace formula or the 
spectral ( function which do not converge in the absolute sense on the real energy 
axis. Therefore, to this date, this approach did not provide an analytic estimate 
of the accuracy of the semiclassical spectrum. 

Under these circumstances, we initiated the present work which addressed 
the problem of the semiclassical accuracy using the approach to be described in 
the sequel. Our main result is that in contrast with the estimate given above, 
the semiclassical error (measured in units of the mean spacing) is independent 
of the dimensionality. Moreover, a conservative estimate of the upper bound 
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for its possible divergence in the semiclassical limit is 0( \ logh\). This is a very 
important conclusion. It allows one to extend many of the results obtained in the 
study of quantum chaos in 2D to higher dimensions, and justifies the use of the 
semiclassical approximation to investigate special features which appear only in 
higher dimensions. We list a few examples of such effects: 

• The dual correspondence between the spectrum of quantum energies and 
the spectrum of actions of periodic orbits [£| [10], [D] was never checked for 
systems in more than 2D. However, if the universality of the quantum spec- 
tral correlations is independent of the number of freedoms, the correspond- 
ing range of correlations in the spectrum of classical actions is expected to 
depend on the dimensionality. Testing the validity of this prediction, which 
is derived by using the trace formula, is of great importance and interest. 
It will be discussed at length in this work. 

• The full range of types of stabilities of classical periodic orbits that includes 
also the loxodromic stability |2j can be manifest only for d > 2. 

• Arnold's diffusion in the KAM regime is possible only for d > 2 (even 
though we do not encountered it in this work). 

Having stated the motivations and background for the present study, we shall 
describe the strategy we chose to address the problem, and the logic behind the 
way we present the results in this report. 

The method we pursued in this first exploration of quantum chaos in 3D, 
was to perform a comprehensive semiclassical analysis of a particular yet typical 
system in 3D, which has a well studied counterpart in 2D. By comparing the 
exact quantum results with the semiclassical theory, we tried to identify possible 
deviations which could be attributed to particular failures of the semiclassical 
approximation in 3D. The observed deviations, and their dependence on h and on 
the dimensionality, were used to assess the semiclassical error and its dependence 
on h. Such an approach requires the assembly of an accurate and complete 
databases for the quantum energies and for the classical periodic orbits. This is 
a very demanding task for chaotic systems in 3D, and it is the main reason why 
such studies were not performed before. 

When we searched for a convenient system for our study, we turned imme- 
diately to billiards. They are natural paradigms in the study of classical and 
quantum chaos. The classical mechanics of billiards is simpler than for systems 
with potentials: The energy dependence can be scaled out, and the system can 
be characterized in terms of purely geometric data. The dynamics of billiards 
reduces to a mapping through the natural Poincare section which is the billiard's 
boundary. Much is known about classical billiards in the mathematical litera- 
ture (e.g. |T2|]), and this information is crucial for the semiclassical application. 
Billiards are also very convenient from the quantal point of view. There are spe- 
cialized methods to quantize them which are considerably simpler than those for 



6 



potential systems |L3| . Some of them are based on the Boundary Integral Method 
(BIM) ||14j| , the KKR method [15|, the scattering approach [17| and various 
improvements thereof |18|, 19, ^Cf| . The classical scaling property is manifest also 
quantum mechanically. While for potential systems the energy levels depend in 
a complicated way on h and the classical actions are non-trivial functions of E, 
in billiards, both the quantum energies and the classical actions scale trivially in 
h and y/E, respectively, which simplifies the analysis considerably. 

The particular billiard we studied is the 3D Sinai billiard. It consists of 
the free space between a 3-torus of side S and an inscribed sphere of radius R, 
where 2R < S. It is the natural extension of the familiar 2D Sinai billiard, and 
it is shown in figure [l] using three complementary representations. The classi- 
cal dynamics consists of specular reflections from the sphere. If the billiard is 
desymmetrized, specular reflections from the symmetry planes exist as well. The 
3D Sinai billiard has several advantages. It is one of the very few systems in 



3D which are rigorously known to be ergodic and mixing pT| , p2| , p3f| . Moreover 



since its introduction by Sinai and his proof of its ergodicity pi], the 2D Sinai 
billiard was subject to thorough classical, quantal and semiclassical investiga- 

Therefore, 



tions 21, 1 



25, 26, 22, 17, 27 



much is known about the 2D 
Sinai billiard and this serves us as an excellent background for the study of the 
3D counterpart. The symmetries of the 3D Sinai billiard greatly facilitate the 
quantal treatment of the billiard. Due to the spherical symmetry of the inscribed 
obstacle and the cubic-lattice symmetry of the billiard (see figure |l](c) ) we are 
able to use the KKR method [^, |2^, |3D|, [T5[ to numerically compute the energy 
levels. This method is superior to the standard methods of computing generic 
billiard's levels. In fact, had we used the standard methods with our present 
computing resources, it would have been possible to obtain only a limited num- 
ber of energy levels with the required precision. The KKR method enabled us 
to compute many thousands of energy levels of the 3D Sinai billiard. The fact 
that the billiard is symmetric means that the Hamiltonian is block-diagonalized 
with respect to the irreducible representations of the symmetry group |JT| . Each 
block is an independent Hamiltonian which corresponds to the desymmetrized 
billiard (see figure |TJ(b) ) for which the boundary conditions are determined by 
the irreducible representations. Hence, with minor changes one is able to com- 
pute a few independent spectra that correspond to the same 3D desymmetrized 
Sinai billiard but with different boundary conditions — thus one can easily ac- 
cumulate data for spectral statistics. On the classical level, the 3D Sinai billiard 
has the great advantage of having a symbolic dynamics. Using the centers of 
spheres which are positioned on the infinite lattice as the building blocks of 
this symbolic dynamics, it is possible to uniquely encode the periodic orbits of 
the billiard [27, 32]. This construction, together with the property that periodic 
orbits are the single minima of the length (action) function [271 |3"2"|, enables us 
to systematically find all of the periodic orbits of the billiard, which is crucial 
for the application of the semiclassical periodic orbit theory. We emphasize that 
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(a) 




(b) (c) 

Figure 1: Three representations of the 3D Sinai billiard: (a) original, (b) 48-fold 
desymmetrized (maximal desymmetrization) into the fundamental domain, (c) 
unfolded to H 3 . 
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performing a systematic search of periodic orbits of a given billiard is far from 
being trivial (e.g. 0, [33], |4], |35|, |36|) and there is no general method of doing so. 



The existence of such a method for the 3D Sinai billiard was a major factor in 
favour of this system. 

The advantages of the 3D Sinai billiard listed above are gained at the expense 
of some problematic features which emerge from the cubic symmetry of the bil- 
liard. In the billiard there exist families of periodic, neutrally stable orbits, the 
so called "bouncing-ball" families that are illustrated in figure |2[. The bouncing- 
ball families are well-known from studies of, e.g., the 2D Sinai and the stadium 



billiards |15| , |17| , p7j , |3q| . These periodic manifolds have zero measure in phase 
space (both in 2D and in 3D), but nevertheless strongly influence the dynamics. 
They are responsible for the long (power-law) tails of some classical distribu- 
tions [[39, 40]. They are also responsible for non-generic effects in the quantum 



spectral statistics, e.g., large saturation values of the number variance in the 2D 



Sinai and stadium billiards j37|. The most dramatic visualization of the effect 
of the bouncing-ball families appears in the function D{1) = ^2 n cos(k n l) — the 
"quantal length spectrum". The lengths I that correspond to the bouncing-ball 
families are characterized by large peaks that overwhelm the generic contribu- 
tions of unstable periodic orbits |38[] (as is exemplified by figure ^). In the 3D 



Sinai billiard the undesirable effects are even more apparent than for the 2D 
billiard. This is because, in general, they occupy 3D volumes rather than 2D 
areas in configuration space and consequently their amplitudes grow as k 1 (to be 
contrasted with k° for unstable periodic orbits). Moreover, for R < S/2 there is 
always an infinite number of families present in the 3D Sinai billiard compared 
to the finite number which exists in the 2D Sinai and the stadium billiards. The 
bouncing balls are thoroughly discussed in the present work, and a large effort 
was invested in devising methods by which their effects could be filtered out. 

After introducing the system to be studied, we shall explain now the way by 
which we present the results. The semiclassical analysis is based on the exact 
quantum spectrum, and on the classical periodic orbits. Hence, the first sections 
are dedicated to the discussion of the exact quantum and classical dynamics in 
the 3D Sinai billiard, and the semiclassical analysis is deferred to the last sections. 
The sections are grouped as follows: 

• Quantum mechanics and spectral statistics (sections § and |3[). 

• Classical periodic orbits (section |J). 

• Semiclassical analysis (sections [|, [| 0). 

In section || we describe the KKR method which was used to numerically 
compute the quantum spectrum. Even though it is a rather technical section, it 
gives a clear idea of the difficulties encountered in the quantization of this system, 
and how we used symmetry considerations and number-theoretical arguments to 
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Figure 2: Some bouncing-ball families in the 3D Sinai billiard. Upper figure: 
Three families parallel to the x, y and z axes. Lower figure: Top view of two 
families. 
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reduce the numerical effort considerably. The desymmetrization of the billiard 
according to the symmetry group is worked out in detail. This section ends with 
a short explanation of the methods used to ensure the completeness and the 
accuracy of the spectrum. 

The study of spectral statistics, section starts with the analysis of the 
integrable billiard (R = 0) case. This spectrum is completely determined by 
the underlying classical bouncing-ball manifolds which are classified according to 
their dimensionality. The two-point form factor in this case is not Poissonian, 
even though the system is integrable. Rather, it reflects the number-theoretical 
degeneracies of the Z 3 lattice resulting in non-generic correlations. Turning to the 
chaotic (R > 0) cases, we investigate some standard statistics (nearest-neighbor, 
number variance) as well as the auto-correlations of the spectral determinant, 
and compare them to the predictions of RMT. The main conclusion of this sec- 
tion is that the spectral fluctuations in the 3D Sinai billiard belong to the same 
universality class as in the 2D analogue. 

Section [| is devoted to the systematic search of the periodic orbits of the 
3D Sinai billiard. We rely heavily on a theorem that guarantees the uniqueness 
of the coding and the variational minimality of the periodic orbit lengths. The 
necessary generalizations for the desymmetrized billiard are also explained. Once 
the algorithm for the computation of periodic orbits is outlined, we turn to the 
definition of the spectrum of lengths of periodic orbits and to the study of its 
statistics. The number of periodic orbits with lengths smaller than L is shown to 
proliferate exponentially. We check also classical sum rules which originate from 
ergodic coverage, and observe appreciable corrections to the leading term due to 
the infinite horizon of the Sinai billiard. Turning our attention to the two-point 
statistics of the classical spectrum, we show that it is not Poissonian. Rather, 
there exist correlations which appear on a scale larger than the nearest spacing. 
This has very important consequences for the semiclassical analysis of the spectral 
statistics. We study these correlations and offer a dynamical explanation for their 
origin. 

The semiclassical analysis of the billiard is the subject of section [5[ As a pre- 
lude, we propose and use a new method to verify the completeness and accuracy 
of the quantal spectrum, which is based on a "universal" feature of the classical 
length spectrum of the 3D Sinai billiard. The main purpose of this section is to 
compare the quantal computations to the semiclassical predictions according to 
the Gutzwiller trace formula, as a first step in our study of its accuracy. Since 
we are interested in the generic unstable periodic orbits rather than the non- 
generic bouncing balls, special effort is made to eliminate the the effects of the 
latter. This is accomplished using a method that consists of taking the derivative 
with respect to a continuous parameterization of the boundary conditions on the 
sphere. 

In section ^ we embark on the task of estimating the semiclassical error of 
energy levels. We first define the measures with which we quantify the semiclas- 
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sical error, and demonstrate some useful statistical connections between them. 
We then show how these measures can be evaluated for a given system using its 
quantal and semiclassical length spectra. We use the databases of the 2D and 3D 
Sinai billiards to derive the estimate of the semiclassical error which was already 
quoted above: The semiclassical error (measured in units of the mean spacing) 
is independent of the dimensionality, and a conservative estimate of the upper 
bound for its possible divergence in the semiclassical limit is 0{\ \ogh\). 

Once we are reassured of the reliability of the trace formula in 3D, we return 
in section [7] to the spectral statistics of the quantized billiard. The semiclassical 
trace formula is interpreted as an expression of the duality between the quantum 
spectrum and the classical spectrum of lengths. We show how the length corre- 
lations in the classical spectrum induce correlations in the quantum spectrum, 
which reproduce rather well the RMT predictions. 

The work is summarized in section |8|. 

To end the introductory notes, a review of the existing literature is in or- 
der. Only very few systems in 3D were studied in the past. We should first 



mention the measurements of 3D acoustic cavities 141, Wa, H3l [44], H5I1 and electro- 



magnetic (microwaves) cavities |47], ^8], f49| . The measured frequency spectra 



were analyzed and for irregular shapes (notably the 3D Sinai billiard) the level 
statistics conformed with the predictions of RMT. Moreover, the length spectra 
showed peaks at the lengths of periodic manifolds, but no further quantitative 
comparison with the semiclassical theory was attempted. However, none of the 
experiments is directly relevant to the quantal (scalar) problem since the acous- 
tic and electromagnetic vector equations cannot be reduced to a scalar equation 
in the configurations chosen. Therefore, these experiments do not constitute a 
direct analogue of quantum chaos in 3D. This is in contrast with flat and thin mi- 
crowave cavities which are equivalent (up to some maximal energy) to 2D quantal 
billiards. 

A few 3D billiards were discussed theoretically in the context of quantum 
chaos. Polyhedral billiards in the 3D hyperbolic space with constant negative 
curvature were investigated by Aurich and Marklof |5(J. The trace formula in 



this case is exact rather than semiclassical, and thus the issue of the semiclassical 
accuracy is not relevant. Moreover, the tetrahedral that was treated had expo- 
nentially growing multiplicities of lengths of classical periodic orbits, and hence 
cannot be considered as generic. Prosen considered a 3D billiard with smooth 
boundaries and 48-fold symmetry [pi], |2(| whose classical motion was almost com- 
pletely (but not fully) chaotic. He computed many levels and found that level 
statistics reproduce the RMT predictions with some deviations. He also found 
agreement with Weyl's law (smooth density of states) and identified peaks of 
the length spectrum with lengths of periodic orbits. The majority of high-lying 
eigenstates were found to be uniformly extended over the energy shell, with no- 
table exceptions that were "scarred" either on a classical periodic orbit or on a 
symmetry plane. Henseler, Wirzba and Guhr treated the iV-sphere scattering 
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systems in 3D [O] in which the quantum mechanical resonances were compared 
to the predictions according to the Gutzwiller trace formula. A good agreement 
was observed for the uppermost band of resonances and no agreement for other 
bands which are dominated by diffraction effects. Unfortunately, conclusive re- 
sults were given only for non-generic configurations of two and three spheres for 
which all the periodic orbits are planar. In addition, it is not clear whether one 
can infer from the accuracy of complex scattering resonances to the accuracy of 
real energy levels in bound systems. Recently, Sieber [52] calculated the 4x4 sta- 
bility (monodromy) matrices and the Maslov indices for general 3D billiards and 
gave a practical method to compute them, which extended our previous results 



for the 3D Sinai billiard 53, p4l 
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2 Quantization of the 3D Sinai billiard 



In the present section we describe the KKR determinant method [28|, |29], |55], [30 



to compute the energy spectrum of the 3D Sinai billiard, and the results of the 
numerical computations. The KKR method, which was used by Berry for the 2D 
Sinai billiard case [ 15| , is most suitable for our purpose since it allows to exploit 



the symmetries of the billiard to reduce the numerical effort considerably. The 
essence of the method is to convert the Schrodinger equation and the boundary 
conditions into a single integral equation. The spectrum is then the set of real 
wavenumbers k n where the corresponding secular determinant vanishes. As a 
matter of fact, we believe that only with the KKR method could we obtain a 
sufficiently accurate and extended spectrum for the quantum 3D Sinai billiard. 
We present in this section also some numerical aspects and verify the accuracy 
and completeness of the computed levels. 

We go into the technical details of the quantal computation because we wish 
to show the high reduction factor which is gained by the KKR method. Without 
this significant reduction the numerical computation would have resulted in only 
a very limited number of levels fl4*6| , 4"8|| . The reader who is not interested in 
these technical details should proceed to subsection pT4| . To avoid ambiguities, 
we strictly adhere to the conventions of 



2.1 The KKR determinant 

We first consider the 3D "Sinai torus" , which is the free space outside of a sphere 
of radius R embedded in a 3-torus of side length S (see fig. p. The Schrodinger 
equation of an electron of mass m and energy E is reduced to the Helmholtz 

equation: 

V 2 ip + k 2 ip = , k = V2mE/h. (1) 

The boundary conditions on the sphere are taken to be the general linear (self- 
adjoint) conditions: 

k cos a ■ ij) + sin a • dnip = , (2) 

where n is the normal pointing outside the billiard, k is a parameter with dimen- 
sions of k, and a G [0, n/2] is an angle that interpolates between Dirichlet (a = 0) 
and Neumann (a = 7r/2) conditions. These "mixed" boundary conditions will 
be needed in section [5] when dealing with the semiclassical analysis. Applying 
the KKR method, we obtain the following quantization condition (see j51J for a 
derivation and for details): 

det [Ai m>Vm >{k) + kPi(kR; k, a)5 w 5 mm >} = , (3) 

l,V = 0,1,2,..., < m < 1 , < m < I' , 
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where k is the wavenumber under consideration and: 



A 



Im.l'm' 



(k) 



4m 1 - 1 ' 



LM 

L = 0,1,2, 



Clm ,lm,l'm' 



M = 0, 



D LM (k) = (-ik) 



h+(kSp)Yl M (^) + -^=S L0 



Clm ,lm,l'm' 

Pi(kR; K, a) 



pez 3 /{0} 

de / dcf)Y LM (e, 4>)Y? m {e, <j>)Y Vml {e, <p) 

o Jo 

kR cos a ■ ni(kR)—kR sin a ■ n[{kR) 
K-Rcosa ■ ji(kR)— kR sin a ■ j[{kR) 
cot[r]i(kR; k, a)] . 



(4) 
(5) 

(6) 

(7) 
(8) 

In the above ji, ni, are the spherical Bessel, Neumann and Hankel functions, 



respectively |5B|, Yi m are the spherical harmonics [5B with argument in the 
direction of p, and rji are the scattering phase shifts from the sphere, subject to 
the boundary conditions @. 

The physical input to the KKR determinant is distributed in a systematic 
way: The terms Ai m ^ m '{k) contain information only about the structure of the 
underlying Z 3 lattice, and are independent of the radius R of the inscribed sphere. 
Hence they are called the "structure functions" |28|, |3(| . Moreover, they depend 
on a smaller number of "building block" functions Dlm(^) which contain the 
infinite lattice summations. The diagonal term kPi(kR)8u>8 mm ' contains the in- 
formation about the inscribed sphere, and is expressed in terms of the scattering 
phase shifts from the sphere. This elegant structure of the KKR determinant (|3|) 
prevails in more general situations and remains intact even if the ~S} lattice is 
replaced by a more general one, or if the "hard" sphere is replaced by a "soft" 
spherical potential with a finite range ("muffin-tin" potential) p8), |30], j29f . This 
renders the KKR a powerful quantization method. In all these cases the struc- 
ture functions Ai m y m i depend only on the underlying lattice, and the relation (|8|) 
holds with the appropriate scattering matrix. Thus, in principle, the structure 
functions (or rather Dlm) can be tabulated once for a given lattice (e.g. cubic) 
as functions of k, and only Pi need to be re-calculated for every realization of the 
potential (e.g. changing R). This makes the KKR method very attractive also 
for a large class of generalizations of the 3D Sinai billiard. 

The determinant (|3|) is not yet suitable for numerical computations. This is 
because the lattice summations in Dlm are only conditionally convergent and 
have to be resummed in order to give absolutely and rapidly convergent sums. 
This is done using the Ewald summation technique, which is described in appen- 
dices |C]-[E]. The further symmetry reductions of the KKR determinant, which 
are one of the most important advantages of this method, are discussed in the 
following. 
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2.2 Symmetry considerations 



As can be seen from equations and from appendix Q the main compu- 

tational effort involved in computing the KKR determinant is consumed in the 
lattice sums Z)lm(^) which need to be evaluated separately for every k. There- 
fore, it is imperative to use every possible means to economize the computational 
effort invested in calculating these functions. For this purpose, we shall exploit 
the cubic symmetry of the 3D Sinai billiard as well as other relations that dras- 
tically reduce the computational effort. 



2.2.1 Group— theoretical resummations 

For the practical (rapidly convergent) computation, the functions Dlm are de- 
composed into three terms which are given in appendix y (see also appendix 
|D|). Equations (|211|) - (|214|) express D^' M as a sum over the direct cubic lattice, 



whereas, D^' is a sum over the reciprocal cubic lattice, which is also a cubic 
lattice. Thus, both sums can be represented as: 

Df M {k) = k )n M {^) , j = i,2. (9) 

We show in appendix |GJ that lattice sums of this kind can be rewritten as: 

d&W = £ £ run**) , (io) 

p p L[Pp) geO h 



where Oh is the cubic symmetry group |3T[, and p p = (11,12,13) resides in the 
fundamental section < %\ < 12 < 23. The terms l(p p ) are integers which are 
explicitly given in appendix |G]. The inner sums are independent of k, and can 
thus be tabulated once for all. Hence the computation of the k dependent part 
becomes 48 times more efficient (for large, finite lattices) when compared to (|9|) 
due to the restriction of p p to the fundamental section. 

A further reduction can be achieved by a unitary transformation from the 
{Ylm} basis to the more natural basis of the irreducible representations (irreps) 
of O h : 

Yt%(V) = Y. a %,MYLM(tt), (11) 

M 

where 7 G [1, . . . , 10] denotes the irrep under consideration, J counts the number 
of the inequivalent irreps 7 contained in L, and K = 1, . . . ,dim(7) is the row 
index within the irrep. The functions Y^j K are known as the "cubic harmonics" 
[ pwjj . Combining ( |T0D and (0), and using the unitarity of the transformation as 



well as the "great orthonormality theorem" of group theory |JT| we arrive at: 

D ( lkk) = J2 a % D< ikk) (12) 
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Dfj{k), = 48^:^^yif(%). (13) 

Pp 

The superscript (s) denotes the totally symmetric irrep. The constant coefficients 
a ij,M can be taken into the (constant) coefficients CLM,im,i'm' resulting in: 

Ai m , Vm >{k) = Am l ~ l 'J2 rLD Lj(k)C LJMilw (14) 

LJ 

D LJ (k) = D L l ]{k)+D L 2 ]{k)+D^{k)5 L » (15) 

CLJ,lm,l'm' = ^sJ*M^LU,lm,Vm! ■ (16) 

M 

We show in appendix [F] that for large L the number of DLj(k)'s is smaller by 
a factor « 1/48 than the number of DlmW's. This means that the entries of 
the KKR determinant are now computed using a substantially smaller number 
of building blocks for which lattice summations are required. Thus, in total, we 
gain a saving factor of 48 2 = 2304 over the more naive scheme (f§-|6|). 

2.2.2 Number— theoretical resummations 

In the above we grouped together lattice vectors with the same magnitude, using 
the geometrical symmetries of the cubic lattice. One can gain yet another re- 
duction factor in the computational effort by taking advantage of a phenomenon 
which is particular to the cubic lattice and stems from number theory. The lengths 
of lattice vectors in the fundamental sector show an appreciable degeneracy, which 
is not connected with the Oh symmetry. For example, the lattice vectors (5, 6, 7) 
and (1, 3, 10) have the same magnitude, a/TTO, and are not geometrically conju- 
gate by Oh- This number-theoretical degeneracy is both frequent and significant, 
and we use it in the following way. Since the square of the magnitude is an integer 
we can write: 



D ( i ) J (k) = f2f if) (p P =V^;k) 



n=l 



V —Y (s) *(n^ 



(17) 



The inner sums incorporate the number theoretical degeneracies. They are k 
independent, and therfore can be tabulated once for all. 

To show the efficiency of ( |17D let us restrict our lattice summation to p p < p max 
(which we always do in practice). For large p max the number of lattice vectors 
in the fundamental domain is 7rp^ ax /36, and the number of summands in (|i~7|) 
is at most p^ax- Thus, the saving factor is at least vrp max /36. In fact, as shown 
in appendix 0, there are only (asymptotically) (5/6)p^ ax terms in (0), which 
sets the saving factor due to number-theoretical degeneracy to be vrp max /30. In 
practice, p max = (9(100) and this results in a reduction factor of about 10, which 
is very significant. 
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2.2.3 Desymmetrization 



The symmetry of the 3D Sinai torus implies that the wavefunctions can be classi- 
fied according to the irreps of Oh |3ll • Geometrically, each such irrep corresponds 



to specific boundary conditions on the symmetry planes that define the desym- 
metrized 3D Sinai billiard (see figure H). This allows us to "desymmetrize" the 
billiard, that is to restrict ourselves to the fundamental domain with specific 
boundary conditions instead of considering the whole 3-torus. We recall that the 
boundary conditions on the sphere are determined by Pi{k) and are independent 
of the irrep under consideration. For simplicity, we shall restrict ourselves to the 
two simplest irreps which are both one-dimensional: 

7 = a: This is the totally antisymmetric irrep, which corresponds to Dirichlet 
boundary conditions on the planes. 

7 = s: This is the totally symmetric irrep, which corresponds to Neumann 
boundary conditions on the planes. 



The implementation of this desymmetrization is straightforward (see [54| for 
details) and results in a new secular equation: 



det 







where 7 is the chosen irrep and: 

LJ 

r»(7) _ \^ J 1 ) J l '> n fon\ 

U LJ,lj,l'j' — ^ a jj,m a 1 j', m ' L/ LJ,lm,l'm' l ZU J 
mm' 

- 2-^ a sJ,M a ij,m a 1 j' , m '^LM,lm,l'm' ■ 



Mmm' 



The desymmetrization of the problem has a few advantages: 

Computational efficiency: In appendix we show that for large L's the num- 
ber of cubic harmonics Y^j K that belong to a one-dimensional irrep is 1/48 
of the number of the spherical harmonics Ylm- Correspondingly, if we trun- 
cate our secular determinant such that L < L max , then the dimension of 
the new determinant (Jl8j) is only 1/48 of the original one @ for the fully 
symmetric billiard. Indeed, the desymmetrized billiard has only 1/48 of the 
volume of the symmetric one, and hence the density of states is reduced by 
48 (for large k). However, due to the high cost of computing a determinant 
(or performing a Singular Value Decomposition) []5H] the reduction in the 
density of states is over-compensated by the reduction of the matrix size, 
resulting in a saving factor of 48. This is proven in appendix [B], where it 
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is shown in general that levels of desymmetrized billiards are computation- 
ally cheaper than those of billiards which possess symmetries. Applied to 
our case, the computational effort to compute a given number N of energy 
levels of the desymmetrized billiard is 48 times cheaper than computing iV 
levels of the fully symmetric billiard. 

Statistical independence of spectra: For each irrep the spectrum is statisti- 
cally independent of the others, since it corresponds to different boundary 
conditions. Thus, if the fully symmetric billiard is quantized, the result- 
ing spectrum is the union of 10 independent spectra (there are 10 irreps 



of Oh ||31||), and significant features such as level rigidity will be severely 
blurred |[59|| . To observe generic statistical properties and to compare with 
the results of RMT, one should consider each spectrum separately, which 
is equivalent to desymmetrizing the billiard. 

Rigidity: The statistical independence has important practical consequences. 
Spectral rigidity implies that it is unlikely to find levels in close vicinity of 
each other. Moreover, the fluctuations in the spectral counting functions 
are bounded. Both features of rigidity are used in the numerical algorithm 
which computes the spectrum, and is described in more detail in section 



2.3 



To summarize this subsection, we have demonstrated that the high symmetry 
features of the 3D Sinai billiard are naturally incorporated in the KKR method. 
This renders the computation of its spectrum much more efficient than in the case 
of other, less symmetric 3D billiards. Thus, we expect to get many more levels 



than the few tens that can be typically obtained for generic billiards J46], H|. In 
fact, this feature is the key element which brought this project to a successful 
conclusion. We note that other specialized computation methods, which were 
applied to highly symmetric 3D billiards, also resulted in many levels [ 191, pQf • 



This completes the theoretical framework established for the efficient numer- 
ical computation of the energy levels. In the following we discuss the outcome of 
the actual computations. 

2.3 Numerical aspects 

We computed various energy spectra, defined by different combinations of the 
physically important parameters: 

1. The radius R of the inscribed sphere (the side S was always taken to be 1). 

2. The boundary conditions on the sphere: Dirichlet / Neumann / mixed: 
< a < it/2. 
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3. The boundary conditions on the symmetry planes of the cube: Dirichlet / 
Neumann. These boundary conditions correspond to the antisymmetric / 
symmetric irrep of Oh, respectively. Due to the lattice periodicity, Dirichlet 
(Neumann) boundary conditions on the symmetry planes induce Dirichlet 
(Neumann) also on the planes between neighbouring cells. 

The largest spectral stretch that was obtained numerically corresponded to R = 
0.2 and Dirichlet boundary condition everywhere. It consisted of 6697 levels in 
the interval < k < 281.078. We denote this spectrum in the following as the 
"longest spectrum". 

The practical application of ( |i~8"|) brings about many potential sources of di- 
vergence: The KKR matrix is infinite dimensional in principle, and each of the 
elements is given as an infinite sum over the cubic lattice. To regulate the infi- 
nite dimension of the matrix we use a physical guideline, namely, the fact that 
for I > kR the phase shifts decrease very rapidly toward zero, and the matrix 
becomes essentially diagonal. Therefore, a natural cutoff is Z max = kR, which is 
commonly used (e.g. |T^|). In practice one has to go slightly beyond this limit, 
and to allow a few evanescent modes: Z max = kR + Z e van- To find a suitable value 
of Zevan we used the parameters of the longest spectrum and computed the 17 
eigenvalues in the interval 199.5 < k < 200 with Z evan = 0,2,4,6,8,10 (Z max has 
to be odd). We show in figure |] the successive accuracy of the computed eigen- 
values between consecutive values of Z e van- The results clearly indicate a 10-fold 
increase in accuracy with each increase of Z C van by 2. A moderately high accuracy 
of O(10~ 4 ) relative to level spacing requires Z e van = 8 which was the value we used 
in our computations. 

To regulate the infinite lattice summations in D LJ we used successively larger 
subsets of the lattice. The increase was such that at least twice as many lattice 
points were used. Our criterion of convergence was that the maximal absolute 
value of the difference between successive computations of Dlj was smaller than 
a prescribed threshold: 

max | Z^-Z^ 1 1 <e. (22) 

The threshold e = 10~ 6 was found to be satisfactory, and we needed to use a 
sub-lattice with maximal radius of 161. 

The KKR program is essentially a loop over k which sweeps the /c-axis in 
a given interval. At each step the KKR matrix M(k) is computed, and then 
its determinant is evaluated. In principle, eigenvalues are obtained whenever 
the determinant vanishes. In practice, however, the direct evaluation of the 
determinant suffers from two drawbacks: 

• The numerical algorithms that are used to compute det M(k) are frequently 
unstable. Hence, it is impossible to use them beyond some critical k which 
is not very large. 
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Figure 3: Accuracy of eigenvalues as a function of the number of evanescent 
modes Z cva n- The case considered was R = 0.2 and Dirichlet boundary conditions 
everywhere. The figure shows the absolute differences of the eigenvalues between 
two successive values of / CV an, multiplied by the smooth level density. That is, 
"0-2" means d(k n ) \k n {l eY3Jl —2) — A;„(/ CV an=0)| = |A7V n |. We show 17 eigenvalues 
in the interval 199.5 < k < 200. 
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• For moderately large fc's, the absolute values of det M(k) are very small 
numbers that result in computer underflows (in double precision mode), 
even for fc-values which are not eigenvalues. 

• Due to finite precision and rounding errors, det M(k) never really vanishes 
for eigenvalues. 

A superior alternative to the direct calculation of the determinant is to use the 



Singular Value Decomposition (SVD) algorithm [58], which is stable under any 
circumstances. In our case, M is real and symmetric, and the output are the 
"singular values" <7; which are the absolute values of the eigenvalues of M. The 
product of all of the singular values is equal to | det M\, which solves the stability 
problem. To cure the other two problems consider the following "conditioning 



measure" : 



dimM(fc) 

r(k) = ^ai(k). (23) 

i=i 

The use of the logarithm circumvents the underflow problem. Moreover, we 
always expect some of the smallest singular values to reflect the numerical noise, 
and the larger ones to be physically relevant. Near an eigenvalue, however, one 
of the "relevant" singular values must approach zero, resulting in a "dip" in the 
graph of r(k). Hence, by tracking function of k, we locate its dips and 

take as the eigenvalues the k values for which the local minima of r are obtained. 
Frequently one encounters very shallow dips (typically -C 1) which are due to 
numerical noise and should be discarded. 

To ensure the location of all of the eigenvalues in a certain k interval, the 
/c-axis has to be sampled densely. However, oversampling should be avoided 
to save computer resources. In order to choose the sampling interval Ak in a 
reasonable way, we suggest the following. If the system is known to be classically 
chaotic, then we expect the quantal nearest-neighbour distribution to follow the 
prediction of Random Matrix Theory (RMT) || . In particular, for systems with 
time reversal symmetry: 

P(s) « ~s , s < 1 , s = (k n+1 - k n ) d{{k n + k n+1 )/2) (24) 

where d(k) is the smooth density of states. The chance of finding a pair of energy 
levels in the interval [s, s + ds] is P{s) ds. The cumulative probability of finding 
a pair in [0, s] is therefore crudely given by: 

J(s)» / P(s') ds' « js 2 , s«l. (25) 
Jo 4 

A more refined calculation, taking into account all the possible relative configu- 
rations of the pair in the interval [0, s] gives: 

Q(s)^^s 2 , (26) 
o 
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If we trace the fc-axis with steps Ak and find an eigenvalue, then the chance that 
there is another one in the same interval Ak is Q(Akd(k)). If we prescribe our 
tolerance Q to lose eigenvalues, then we should choose: 

m = *1 (27) 

d[k) d[k) V 7r 

In the above, we assumed that the dips in r(k) are wide enough, such that they 
can be detected over a range of several Afc's. If this is not the case and the dips 
are very sharp, we must refine Ak. In our case dips were quite sharp, and in 
practice we needed to take Q of the order 10~ 5 -j- 10 -6 . 



2.4 Verifications of low— lying eigenvalues 

After describing some numerical aspects of the computation, we turn to various 
tests of the integrity and completeness of the computed spectra. In this subsection 
we compare the computed low-lying eigenvalues for R > with those of the R = 
case. In the next one we compare the computed stair-case function to Weyl's 
law. 

The theoretical background for the comparison between low-lying eigenvalues 
to those of the R = case is as follows. The lowest I value, for which there exist 
antisymmetric cubic harmonics, is I = 9 [57||. Consequently, for cases with Dirich- 



let conditions on the symmetry planes, the lowest /-values in the KKR matrix is 
/ = 9. Thus, for kR < 9 the terms Pi(kR) in equation ([18]) are very small, and 
the matrix is essentially as if the inscribed sphere was not present. In that case 
of the "empty tetrahedron" the eigenvalues can be calculated analytically: 



2tt 



k*=° = — V/ 2 + m 2 + n 2 , < / < m < n . (21 
S 



k n ^k^=° for k n <9/R. (29) 



We hence expect: 



Similar considerations were used by Berry [yj| for the 2D Sinai billiard, where 
he also calculated the corrections to the low-lying eigenvalues. In figure [| we 
plot the unfolded difference AN n = d(k n ) \ k n — k^ =0 \ for the longest spectrum 
(R = 0.2, Dirichlet everywhere). One clearly observes that indeed the differences 
are very small up to k = 9/0.2 = 45, and they become of order 1 afterwards, as 
expected. This confirms the accuracy and completeness of the low-lying levels. 
Moreover, it verifies the correctness of the rather complicated computations of 
the terms Ai^vy which are due to the cubic lattice. 
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Figure 4: The unfolded differences AN n for the low-lying levels of the 3D Sinai 
billiard with R = 0.2 and Dirichlet everywhere. We indicated by the vertical line 
k = 45 the theoretical expectation for transition from small to large AiV. The 
line AiV = was slightly shifted upwards for clarity. 
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2.5 Comparing the exact counting function with Weyl's 
law 



It is by now a standard practice (see e.g. |I7|) to verify the completeness of a 
spectrum by comparing the resulting stair-case function N(k) = #{k n < k} 
to its smooth approximation N(k), known as "Weyl's law". In appendix p] we 
derived Weyl's law for the 3D Sinai billiard (equation Q290Q ), and now consider 
the difference N osc (k) = N(k) — N(k). Any jump of N osc by ±1 indicates a 
redundant or missing eigenvalue. In fact, this tool is of great help to locating 
missing eigenvalues. In figure |^ we plot N osc for the longest spectrum. It is 
evident that the curve fluctuates around with no systematic increase/decrease 
trends, which verifies the completeness of the spectrum. The average of N osc 
over the available fc-interval is (—4) • 10~ 4 which is remarkably smaller than any 
single contribution to N (note that we had no parameters to fit). This is a very 
convincing verification both of the completeness of the spectrum as well as the 
accuracy of the Weyl's law ( |290| ). We also note that the typical fluctuations grow 
quite strongly with k. This is due to the effects of the bouncing-ball families (see 
section p]) and will be discussed further in section jO . 
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Figure 5: N osc (k) for the longest spectrum of the 3D billiard. The data are 
smoothed over 50 level intervals. 
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3 Quantal spectral statistics 



Weyl's law predicts the smooth behaviour of the quantal density of states. There 
is a wealth of information also in the fluctuations, and their investigation is 
usually referred to as "spectral statistics". Results of spectral statistics that 
comply with the predictions of Random Matrix Theory (RMT) are generally 
considered as a hallmark of the underlying classical chaos |24], |59], |2|, |60], [17j . 

In the case of the Sinai billiard we are plagued with the existence of the non- 
generic bouncing-ball manifolds. They influence the spectral statistics of the 3D 
Sinai billiard. It is therefore desirable to study the bouncing balls in some detail. 
This is done in the first subsection, where we discuss the integrable case (R = 0) 
that contains only bouncing-ball manifolds. 

For the chaotic cases R > we consider the two simplest spectral statistics, 
namely, the nearest-neighbour distribution and two-point correlations. We com- 
pute these statistics for the levels of the 3D Sinai billiard, and compare them 
to RMT predictions. In addition, we discuss the two-point statistics of spectral 
determinants that was recently suggested by Kettemann, Klakow and Smilansky 
|61|j as a characterization of quantum chaos. 



3.1 The integrable R = case 

If the radius of the inscribed sphere is set to 0, we obtain an integrable billiard 
which is the irreducible domain whose volume is 1/48 of the cube. It is plotted 
in figure |^. This tetrahedron billiard is a convenient starting point for analyzing 
the bouncing-ball families, since it contains no unstable periodic orbits but only 
bouncing balls. Quantum mechanically, the eigenvalues of the tetrahedron are 
given explicitly as: 

2?r 



k(nmi) = —Vn 2 + m 2 + P , 0<n<m</GlN. (30) 

The spectral density dR =0 (k) = J2o<n<m<i $(.k — k( nm [)) can be Poisson resummed 
to get: 

S 3 k 2 



d R=0 (k) = ^ sine (kS a/ 'p 2 + q 2 + r 2 ) 



_l_ J^iL\^ cos (kSp) -\ cos | k-^=p 




+ ^E cos (^) • ( 31 ) 



s 

6^ 
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In the above sinc(x) = sin(x)/x, sinc(O) = 1, and Jo is the zeroth order Bessel 
function. Let us analyze this expression in some detail. Terms which have all 
summation indices equal to give the smooth part of the density, and all the 
remaining terms constitute the oscillatory part. Collecting the smooth terms 
together we get: 

e3j,2 c2j, q c; 

= ^ - 3^(1 + V2) + ^(27 + 9V2 + SVS) - -5(k - 0) . (32) 



This is Weyl's law for the tetrahedron, which exactly corresponds to ( |290| ) with 
R = (except the last term for which the limit — ^ is different). 

As for the oscillatory terms, it is first useful to replace Jq{x) by its asymptotic 
approximation which is justified in the semiclassical limit k — ► oo: 

Jo(x) ~ \ — cos (x — — ) , x — > oo . (33) 



TCX \ 4- 

Using this approximation we observe that all of the oscillatory terms have phases 
which are of the form (k x length + phase). This is the standard form of a 
semiclassical expression for the density of states of a billiard. To go a step further 
we notice that the leading-order terms, which are proportional to k l (first line of 
(|3"TD), have lengths of S\Jp 2 + q 2 + r 2 which are the lengths of the periodic orbits 
of the 3-torus, and therefore of its desymmetrization into the tetrahedron. This 
conforms with the expressions derived by Berry and Tabor ||| |64| for integrable 



systems. The other, sub-leading, oscillatory contributions to fl3T|) correspond to 
"improper" periodic manifolds, in the sense that their dynamics involves non- 
trivial limits. Some of these periodic orbits are restricted to symmetry plane or 
go along the edges. Of special interest are the periodic orbits that are shown in 
figure They are isolated, but are neutrally stable and hence are non-generic. 
Their contributions are contained in the last two terms of (|31|), and the one with 
length S/y/3 is the shortest neutral periodic orbit. Other sub-leading oscillatory 
contributions are discussed in p4| . We therefore established an interpretation 
in terms of (proper or improper) classical periodic orbits of the various terms of 



3.1.1 Two-point statistics of the integrable case 

We continue by investigating the two-point statistics of the tetrahedron, which 
will be shown to provide some nontrivial and interesting results. Since we are 
interested in the limiting statistics as k —>■ oo we shall consider only the leading 
term of (|3T|), which is the first term. Up to a factor of 48, this is exactly the 
density of states d T z of the cubic 3-torus, and thus for simplicity we shall dwell 
on the 3-torus rather than on the tetrahedron: 

d T3 (k) =£^-fp)=^£ smc(kSp) . (34) 
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Figure 6: Upper: Geometry of the tetrahedron (R = 0) billiard. Lower: Neutral 
periodic orbits in the desymmetrized 3D Sinai. The billiard is indicated by bold- 
face edges. Dot-dash line: The shortest neutral periodic orbit of length S/ a/3. 
Double dot-dash line: Neutral periodic orbit of length S/ \/2. 
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We observe that both the quantal spectrum and the classical spectrum (the set of 
lengths of periodic orbits) are supported on the cubic lattice Z 3 , and this strong 
duality will be used below. 

The object of our study is the spectral form factor, which is the Fourier 



transform of the two-point correlation function of the energy levels |p9| . For 
billiards it is more convenient to work with the eigenwavenumbers k n rather than 
with the eigenenergies E n . Here the form factor is given by: 



K(r; k) = 1 



7! 2 



exp [27iid(k)k n r~j 



n=ni 



(35) 



In the above N = rii — ri\ + 1, and k n are the eigenvalues in the interval [k ni , k n2 ] 
centered around k = (k ni + k n2 )/2. It is understood that the interval contains 
many levels but is small enough such that the average density is almost a constant 
and is well approximated by d(k). 

In the limit r — ► oo the phases in (|35|) become random in the generic case, 
and therefore K(t) — > 1. However, if the levels are degenerate, more care should 
be exercised, and one obtains: 

K ( r ; k ) = jf}29k{k n ) = -j^Yl ' T ~* 00 ' ( 36 ) 

n i 

where gk{k n ) is the degeneracy of k n and the primed sum is only over distinct 
values of fcj. Since N = Yl'i9k(h) we obtain: 

/C(r;*) = S^M = (£|W>, T ^ 00 , (37 ) 
J2i9k{ki) {9k{k)) 

where (■) denotes an averaging over fc^'s near k. In the case of a constant g the 
above expression reduces to K(t — > oo) = g, but it is important to note that 
K{t — ► oo) 7^ (g) for non-constant degeneracies. Using the relation p = kS/ (2tt) 
(see equation (|3^) ) and equations ( p74j ), (|276|) in appendix [H] we get: 



(g}(kS/(2ir))) PS, 

Mr;fc) = pW = ^ fc ' r ^°°' (38) 

where /3 ~ 9.8264 is a constant. That is, contrary to the generic case, the satu- 
ration value of the form factor grows linearly with k due to number-theoretical 
degeneracies. 

Turning to the form factor in the limit r — > 0, we first rewrite fl3~4]) as dx3(k) = 
d(k) + J2j Aj sm(kLj). Then, using the diagonal approximation as suggested by 
Berry [|J p3J , and taking into account the degeneracies ge(Lj) of the lengths we 
have: 

K{r; k) = 4d^)^ 2(Lj)|Ai|25(T " Lj/Lh) ' T ^ 1 ■ (39) 
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In the above the prime denotes summation only over distinct classical lengths, 
and Lh = 2ird(k) is called the Heisenberg length. The coefficients Aj are func- 
tions of Lj and therefore can be replaced by the function A(t). For r large 
enough such that the periodic manifolds have a well-defined classical density 
d c \(£), the summation over delta functions can be replaced by multiplication with 
-^H^ci(^) / (ge{fy) with I = Lht such that: 

A straightforward calculation shows that the term in brackets is simply 1, which 
is the generic situation for the integrable case (Poisson statistics) ||, |66j. Hence, 
we obtain: 

* (T:i) = S- t -°- (41) 

Since, as we noted above, the lengths of the classical periodic orbits are supported 
on the Z 3 lattice, we can write using £ = Sp: 



2 



* ( ^)=<M=£*-r, r^O. (42) 



WIS)) 



7T 



where we used again equations fl274j ), (|276|) . This is a very surprising result, since 



it implies that contrary to the generic integrable systems, which display Poisson 
level statistics with K — 1, here K oc r which is typical to chaotic systems! This 
peculiarity is manifestly due to the number theoretical degeneracies of Z 3 . 

If we now combine the two limiting behaviours of the form factor in the 
simplest way, we can express it as a scaled RMT-GUE form factor: 

K T3 (r;k) « • # gue ( 7 t) (43) 

where = S/3k/(2ir) and 7 = 2Sk. For the tetrahedron we have the same 
result with — ► Koc/48 and 7 — > 7/48. This prediction is checked and ver- 
ified numerically in figure [7] where we computed the quantal form factor of the 
tetrahedron around various fc-values. The agreement of the two asymptotes to 
the theoretical prediction (|43D is evident and the difference from Poisson is well 
beyond the numerical fluctuations. 



3.2 Nearest— neighbour spacing distribution 

We now turn to the chaotic case R > 0. One of the most common statistical 
measures of a quantum spectrum is the nearest-neighbour distribution P(s). 
If fact, it is the simplest statistics to compute from the numerical data. We 
need only to consider the distribution of the scaled (unfolded) spacings between 
neighbouring levels: 

s n = N(k n+1 ) - N(k n ) w d(k n )(k n+1 - k n ) . (44) 
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Figure 7: The scaled quantal form factor of the tetrahedron for various /c-values 
compared with GUE and Poisson. Note the log-log scales. 
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It is customary to plot a histogram of P(s), but it requires an arbitrary choice of 
the bin size. To avoid this arbitrariness, we consider the cumulant distribution: 



J(s) = / ds'P(s') (45) 



for which no bins are needed. Usually, the numerical data are compared not 
to the exact -Prmt(s) but to Wigner's surmise 0, which provides an accurate 
approximation to the exact Prmt(s) in a simple closed form. In our case, since we 
found a general agreement between the numerical data and Wigner's surmise, we 
choose to present the differences from the exact expression for Jgoe( s ) taken from 
Dietz and Haake ||67|| . In figure |8] we show these differences for R = 0.2, 0.3 and 
Dirichlet boundary conditions (6697 and 1994 levels, respectively). The overall 
result is an agreement between the numerical data and RMT to better than 
4%. This is consistent with the general wisdom for classically chaotic systems in 



lower dimensions, and thus shows the robustness of the RMT conjecture |24j for 
higher-dimensional systems (3D in our case). 

Beyond this general good agreement it is interesting to notice that the dif- 
ferences between the data and the exact GOE for R — 0.2 seem to indicate 
a systematic modulation rather than a statistical fluctuation about the value 
zero. The same qualitative result is obtained for other boundary conditions with 
R = 0.2, substantiating the conjecture that the deviations are systematic and 
not random. For R = 0.3 the differences look random and show no particular 
pattern. However, for the upper third of the spectrum one observes structures 
which are similar to the R = 0.2 case (see figure f| lower part). 

Currently, we have no theoretical explanation of the above mentioned system- 
atic deviations. They might be due to the non-generic bouncing balls. To assess 
this conjecture we computed P(s) for R = 0.2, 0.3 with Dirichlet boundary con- 
ditions in the spectral interval 150 < k < 200. The results (not shown) indicate 
that the deviations are smaller for the larger radius. This is consistent with the 
expected weakening of the bouncing-ball contributions as the radius grows, due 
to larger shadowing and smaller volumes occupied by the bouncing-ball families. 
Hence, we can conclude that the bouncing balls are indeed prime candidates for 
causing the systematic deviations of P(s). It is worth mentioning that a detailed 
analysis of the P{s) of spectra of quantum graphs show similar deviations from 
Prmt(s) ||. 



3.3 Two— point correlations 

Two-point statistics also play a major role in quantum chaos. This is mainly due 
to their analytical accessibility through the Gutzwiller trace formula as demon- 
strated by Berry [HI |5"5|. There is a variety of two-point statistical measures 



which are all related to the pair-correlation function [59]. We chose to focus on 



E 2 (/) which is the local variance of the number of levels in an energy interval 
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Figure 8: Differences of integrated nearest-neighbour distribution for R = 0.2 
(up) and R = 0.3 (down). Set #1, 2, 3 refer to the division of the spectrum into 
3 domains. Data are slightly smoothed for clarity. 
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that has the size of / mean spacings. The general expectation for generic sys- 
tems, according to the theory of Berry ||, |65|, is that £ 2 should comply with 
the predictions of RMT for small values of I (universal regime) and saturate to 
a non-universal value for large Z's due to the semiclassical contributions of short 
periodic orbits. The saturation value in the case of generic billiards is purely 
classical (fc-independent). The effect of the non-generic bouncing-ball manifolds 
on two-point spectral statistics was discussed in the context of 2D billiards by 
Sieber et al. |57| (for the case of the stadium billiard). They found that £ 2 can 



be decomposed into two parts: A generic contribution due to unstable periodic 
orbits and a non-generic contribution due to bouncing balls: 

S 2 (0«E 2 UPO (0 + ^ b (0- (46) 

The term E bb has the structure: 

£ 2 b (Z) = kF stsdium (l/d(k)), (47) 

where -Fstadium is a function which is determined by the bouncing balls of the 
stadium billiard, and is given explicitly in ||37| . In particular, for large values of 



I the term E bb fluctuated around an asymptotic value: 

^bb(0 ~ ^ -^stadium 

(oo) , / — > OO . (48) 

One can apply the arguments of Sieber et al. [37] to the case of the 3D Sinai 



billiard and obtain for the leading order bouncing balls (see ([34])): 

-^3Dsb (WO) , (49) 

with F 3 Dsb characteristic to the 3D Sinai billiard. Asymptotically, we expect: 

£ 2 b (Z) w fc 2 F 3Dsb (oo) , l^oo. (50) 

The function -F 3 D S b can be written down, albeit it contains the areas of the cross- 
sections of the various bouncing-ball manifolds, for which we have no explicit 
expressions. Therefore, we shall investigate the scaling features of S bb without 
insisting on its explicit from. 

The numerical computations of S 2 for the longest spectrum (R = 0.2, Dirich- 
let everywhere) are shown in figure |9|. We divided the spectrum into 4 intervals 
such that d did not vary much within each interval. This is a pre-requisite for 
a meaningful semiclassical analysis. It is evident from the figure that for small 
values of I (up to ~ 1) there is an agreement with GOE. Moreover, the agreement 
with GOE is much better than with either GUE or Poisson, as expected. This is 
in agreement with the common knowledge in quantum chaos and again, sub- 



stantiates the RMT conjecture also for chaotic systems in 3D. For larger / values 
there are marked deviations which saturate into oscillations around a /c-dependent 
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asymptotic values. It is clearly seen that the saturation values grow faster than 
k, which is consistent with (0). To test ( f49|) quantitatively, we plotted in figure 



10 the rescaled function: 



S L(i; k) = ^ [X 2 (qd(k)) - E 2 GOE (qd(k))] (51) 

which according to ((49]) is the fc-independent function -F 3 Dsb(<?)- Indeed, there is 
a clear data collapse for q < 5, and the saturation values of Sl h are of the same 



magnitude for all values of k. This verifies (J49 ) and demonstrates the important 



part which is played by the bouncing balls in the two-point (long range) statistics. 

For generic systems the agreement between S 2 and RMT should prevail up 
to I*, where: 

L H (k) 2ird(k) 

1 = —r = —j • (52) 

In the above Lh is Heisenberg length and L min is the length of the shortest 
periodic orbit. For the cases shown in figure |9] the value of I* is of the order 
of 100. Nevertheless, the deviations from the universal predictions start much 
earlier. This is again a clear sign of the strong effect of the bouncing-balls. To 
substantiate this claim, we compare in figure [ll| the number variances for R = 0.2 
and R = 0.3 in the same k interval and with the same boundary conditions 
(Dirichlet). The influence of the bouncing-balls is expected to be less dominant 
in the R = 0.3 case, since there are fewer of them with smaller cross sections. 
This is indeed verified in the figure: The agreement with GOE predictions lasts 
much longer (up to I ~ 6) in the R = 0.3 case, and the saturation value is smaller, 
as expected. 

3.4 Auto-correlations of spectral determinants 

The two-point correlations discussed above are based on the quantal spectral den- 
sities. Kettemann, Klakow and Smilansky |61] introduced the auto-correlations 



of quantal spectral determinants as a tool for the characterization of quantum 
chaos. Spectral determinants are defined as: 

Z(E) =0 E = E n , (53) 

that is, they are iff E is an eigenenergy. The (unnormalized) correlation function 
of a spectral determinant is defined as: 



1 pE+AE/2 

c(Wi E) s ~ke Je—ae/2 &ez ( e + a) z " i E ' - a) • w « AE ■ (54) 



rE+AE/2 
'E-AE/2 

There are various motivations to study the function C{u) []6l| : 



37 



0.0 



120 < k < 160 
160 < k < 200 
200 < k < 240 
240 < k < 280 



10 



15 



20 



Figure 10: Rescaled number variance ([y]) for the longest spectrum. 
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Figure 11: Comparison between the number variances for two different radii 
R = 0.2, 0.3 of the inscribed sphere of the 3D Sinai billiard. In both cases we 
considered the spectral interval 120 < k < 160 and used Dirichlet boundary- 
conditions. 



39 



1. There is a marked difference in the behaviour of C(u>) for rigid and non- 
rigid spectra. For completely rigid spectra the function C(u) is oscillatory, 
while for Poissonian spectra it rapidly decays as a Gaussian. For the RMT 
ensembles it shows damped oscillations which are due to rigidity. 

2. The function C(u) contains information about all n-point correlations of 
the spectral densities. Thus, it is qualitatively distinct from the two-point 
correlations of spectral densities and contains new information. 

3. The Fourier transform of C(u) exhibit in an explicit and simple way sym- 
metry properties which are due to the reality of the energy levels. 

4. In contrast to spectral densities, the semiclassical expressions for spectral 
determinants can be regularized using the method of Berry and Keating 



59 1 . Regularized semiclassical spectral determinants contain a finite num- 



ber of terms, and are manifestly real for real energies. 

5. The semiclassical expression for C(u) is closely related to the classical Ru- 
elle zeta function. 

To study C{uS) numerically, regularizations are needed. For the 3D Sinai 
billiard the longest spectrum was divided into an ensemble of 167 intervals of 
N = 40 levels, and each interval was unfolded to have mean spacing 1 and was 
centered around E = 0. For each unfolded interval Ij the function Cj(u) was 



computed using equation (69) of plj , with AE = yN. The ensemble average 
function C(u) was normalized such that C(0) = 1. The results of the computation 



are shown in figure 12. The agreement with RMT is quite good up to uo 



that is for short energy scales for which we indeed expect universality to hold. 
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Figure 12: The two-point correlation function of spectral determinants C(u) for 
the 3D Sinai billiard (longest spectrum). The spectrum was divided into 167 
intervals of 40 levels each and the average correlation function is shown. The 
continuous line is the RMT-GOE theoretical curve, and the dashed line is the 
numerical correlation. The correlation function is normalized to 1 for uj = 0. 



With kind permission from the authors of [61 



41 



4 Classical periodic orbits 



In this section we present a comprehensive study of the periodic orbits of the 
3D Sinai billiard. By "periodic orbits" we mean throughout this section generic, 
isolated and unstable periodic orbits which involve at least one bounce from the 
sphere. Thus, bouncing-ball orbits are not treated in this section. The classical 
periodic orbits are the building blocks for the semiclassical Gutzwiller trace for- 
mula, and are therefore needed for the semiclassical analysis to be presented in 
the next section. 

4.1 Periodic orbits of the 3D Sinai torus 

We found it necessary and convenient to first identify the periodic orbits of the 
symmetric 3D Sinai billiard on the torus, and to compute their lengths and 
stabilities. The periodic orbits of the desymmetrized 3D Sinai billiard could 
then be derived by an appropriate classical desymmetrization procedure. 

The basic problem is how to find in a systematic (and efficient) way all the 
periodic orbits of the 3D Sinai billiard up to a given length L max . In dealing 
with periodic orbits of the Sinai billiard it is very helpful to consider its unfolded 
representation that tessellates 1R 3 — as is shown in figure [I| We start by con- 
sidering the periodic orbits of the fully symmetric 3D Sinai billiard on the torus 
(ST). This case is simpler than the desymmetrized billiard, since it contains no 
boundaries and the tiling of the H 3 space is achieved by simple translations along 
the cubic lattice Z 3 . In the unfolded representation every orbit is described by a 
collection of straight segments which connect spheres. At a sphere, the incident 
segment reflects specularly. A periodic orbit of period n is not necessarily peri- 
odic in the unfolded representation, but rather, it obeys the restriction that the 
segments repeat themselves after n steps modulo a translation by a lattice vector 
(see figure |T3"D . If we fix an origin for the lattice, we can assign to every orbit (not 
necessarily periodic) a "code word" by concatenating the "addresses" (locations 
of the centers on the Z 3 lattice) of the spheres from which it reflects. The code 
word can consist of either the absolute addresses of the spheres or alternatively, 
the address of the sphere relative to the previous one. We shall adopt the latter 
convention and use the relative addresses as the "letters" from which the code 
word is composed. This relative coding has the advantage that a periodic orbit 
is represented by a periodic code word. The number of possible letters ( "alpha- 
bet") is obviously infinite and the letter (0, 0, 0) is excluded. A periodic orbit can 
be represented by any cyclic permutation of its code. To lift this ambiguity, we 
choose a convenient (but otherwise arbitrary) lexical ordering of the letters and 
use the code word which is lexically maximal as the unique representative of the 
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periodic orbit: 

(periodic orbit of ST) i — > W = (wi, W2, ■ ■ ■ , w n ) , Wi G Z 3 \(0, 0, 0) 
W = max{W } PW, P 2 W, . . . , P n ~ l W} , 



(55) 



where PW = (u>2, W3, . . . ,w n ,wi) is the operation of a cyclic permutation of the 
code word. 

Let us consider the code word W with n letters: 

W = {w 1 ,w 2 , . . . ,w n ) , Wi = (w ix ,w iy ,w iz ) . (56) 

It relates to the n + 1 spheres centered at c\ = (0,0,0), C2 = u>i, C3 = u>i + 
u>2, • • • , c n+ i = lOi + • • • + w n . Let us choose arbitrary points on each of the 
spheres, and connect them by straight segments. We get a piecewise straight 
line which leads from the first to the last sphere, which, in general, is not a 
classical orbit because the specular reflection conditions are not satisfied. To find 
a periodic orbit, we specify the positions of the points on each sphere by two 
angles 9{ , The length of the line is a function of {(8i, </?i)|i=i,-,n}- Periodic 
orbits on the ST must have identical coordinates for the first and the last points 
(modulo a lattice translation), hence 9 n+ i = 81, f n+ i = <pi and we have only 
2n independent variables to completely specify a periodic set of segments, with 
length: 



where Lj are the lengths of the segments that correspond to the letter W{. To sat- 
isfy the condition of specular reflection we require that the length Lw is extremal 
with respect to any variation of its variables. 

The following theorem guarantees two essential properties of the coding and 
of the periodic orbits which are identified as the extrema of ([57]) |27|, p2| : 



Theorem: To each code word W of the 3D ST there corresponds at most one 
periodic orbit which is the only minimum of Ly/. 

The theorem contains two statements: First, that periodic orbits are necessarily 
minima of the length, and not saddles or maxima. Second, that there are no 
local minima besides the global one. The phrase "at most" in the theorem above 
needs clarification: For each code word W the length function Lw is a continuous 
function in all of its variables over the compact domain which is the union of the 
spheres. Therefore Lw must have a global minimum within this domain. This 
minimum can be, however, classically forbidden, meaning that at least one of its 
segments cuts through one or more spheres in the lattice (that might or might 
not be a part of the code) rather than reflecting from the outside. This is called 



"shadowing". An example is shown in figure IT3. The forbidden periodic orbits 
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are excluded from the set of classical periodic orbits. (They also do not contribute 



to the leading order of the trace formula JTTJ, [L5| and therefore are of no interest 
in our semiclassical analysis.) If all the segments are classically allowed, then 
we have a valid classical periodic orbit. Finally we would like to mention that 



the minimality property was already implied in the work of Sieber |60|| , and 
the explicit versions of the theorem were proved simultaneously by Bunimovich 



(general formulation, applies in particular to the 3D case) and Schanz |p2 
(restricted to the 2D Sinai billiard). 

The number of letters in the codes of periodic orbits of length less than L max 
can be bounded from above by the following argument. To each letter w there 
corresponds a minimal segment length L min (w) > which is the minimum dis- 
tance between the spheres centered at (0,0,0) and at w = (w x ,w y ,w z ): 



L min (w) = S^Jwl + w 2 y + wl-2R. {hi 



In the above, S is the lattice constant (torus's side) and R is the radius of the 
sphere. The smallest possible L min (w) is obtained for w = (1,0,0) and equals 
S — 2R = L min . We readily conclude that the code word cannot contain more 
letters than the integer part of L max / L min . 

We are now in a position to formulate an algorithm for a systematic search 
of all the periodic orbits of length up to -L max of the 3D Sinai torus: 

1. Collect all of the admissible letters into an alphabet. An admissible letter 
w satisfies: 

(a) w + (0,0,0). 

(b) w is not trivially impossible due to complete shadowing, e.g., like 
(2,0,0) =2 x (1,0,0). 

(c) L m [ n (w) < -^max- 

2. Define an arbitrary lexical order of the letters. 

3. From the admissible alphabet construct the set of admissible code words 
W = (u>i, . . . , w n ), such that: 

(a) L min (W) = J27=l L mm( w i) < ^rnax- 

(b) Wi 7^ w i+ i — no a-priori complete shadowing. 

(c) W is lexically maximal with respect to cyclic permutations: W = 
maxjPW, i — 0, . . . , n — 1}. 

4. For each candidate code word W minimize numerically the function Lw- 
According to the theorem, there should be exactly one minimum, which is 
the global one. 
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5. Check whether the resulting periodic orbit is shaded. Accept only periodic 
orbits which are not shaded. 

Once the periodic orbit is identified, its monodromy (stability) matrix is com- 
puted according to the recipe given in appendix [J]. 

4.2 Periodic Orbits of the 3D Sinai billiard — Classical 
desymmetrization 

If we desymmetrize the ST into the Sinai billiard (SB), we still find that the 
SB tessellates the H 3 space. Hence, each periodic orbit of the ST is necessarily 
also a periodic orbit of the SB. The converse is not true, i.e., periodic orbits of 
the SB are not necessarily periodic in ST. However, it is easy to be convinced 
that if a periodic orbit of SB is repeated sufficiently many times, it becomes also 
periodic in ST. An example is shown in figure pl| From a more abstract point 
of view, this is because the cubic group Oh is finite. Thus in principle one could 
use the algorithm given above to systematically find all the periodic orbits of 
the SB. This is, however, highly inefficient because by analyzing the group Oh 
we find that in order to find all the periodic orbits of the SB up to L max we 
must find all of the periodic orbits of ST up to 6L max . Due to the exponential 
proliferation of periodic orbits this would be a colossal waste of resources which 
would diminish our ability to compute periodic orbits almost completely. To 
circumvent this difficulty, without losing the useful uniqueness and minimality 
properties which apply to the ST, we make use of the property that periodicity 
in the SB is synonymous to periodicity in ST modulo an element g G Oh- This 
simple geometrical observation is a manifestation of the fact that the tiling of 
H 3 by the SB is generated by the group Oh <S> 2 3 . Thus, we can represent the 
periodic orbits of the SB by using their unfolded representation, augmented by 
the symmetry element g according to which the periodic orbits closes: 

Periodic orbit of SB i — > W = (W; g) = (wi,w 2 , ■ ■ ■ , w n ; g) . (59) 

The coding is not yet well-defined since a given periodic orbit can in general 
be represented by several codes. Similarly to the case of the ST, there is a 
degeneracy with respect to the starting point. However, in the case of the SB 
this is not simply related to cyclic permutations. Rather, if a periodic orbit is 
described by (wi,w 2 , . . . ,w n ;g) then it is also described by: 

(w 2 , w 3 , . . . , w n , gwi, g), (w 3 , tu 4 , • • • , gwi, gw 2 ; g),..., 
(gw u gw 2 , . . . , gw n ; g), (gw 2 , gw 3l . . . , g 2 wf, </),... 

{$<m-i Wu gm-i W2j . . . 5 g^)-^ Wn . . . . 5 ($m-i Wni Wl ,w 2 ,..., w n _ i; g) . 

(60) 
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Figure 14: A shadowed (classically forbidden) periodic orbit of the Sinai 3-torus. 




Figure 15: Desymmetrization of orbits from the Sinai torus to the Sinai billiard. 
For clarity we show an example in 2D. Left: A primitive periodic orbit in the ST. 
Right: The corresponding periodic orbit in the SB. We observe that the latter is 
4 times shorter than the former. 
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In the above (f>(g) is the period of g, which is defined as the smallest natural 
number for which g^^ = e, where e is the identity operation. For Oh in particular 
<f)(g) G {1,2,3,4,6}. The above generalized cyclic permutation invariance is due 
to the periodicity modulo g of the periodic orbits of the SB in the unfolded 
representation. In addition to the generalized cyclic invariance there is also a 
geometrical invariance of orbits of the SB in the unfolded representation. Indeed, 
if we operate on an orbit in the unfolded representation with any h G Oh we 
obtain the same orbit in the SB. This symmetry is carried over also to the codes. 
If a periodic orbit is described by (wi, w 2 , • • • , w n ; g) then it is also described by: 

(hwi,hw 2 ,...,hw n ;hgh~ 1 ) MhEOh- (61) 

To summarize, a periodic orbit of the SB can be encoded into a code word up to 
degeneracies due to generalized cyclic permutations and geometrical operations. 
The set of operations which relate the various codes for a given periodic orbit is 
a group to which we refer as the invariance group. 

In order to lift this degeneracy and to obtain a unique mapping of periodic 
orbits of the SB to code words we need to specify a criterion for choosing exactly 
one representative. There are many ways of doing this, but we found it convenient 
to apply the natural mapping of periodic orbits of the SB to those of the ST, and 
there, to choose the maximal code. More specifically: 

1. Select the alphabet according to the rules prescribed in the preceding sub- 
section, and define ordering of letters. 

2. Extend the word W into W: 

W = (w 1 ,w 2 ,...,w n ,gw 1 ,gw 2 ,...,gw n , 

g 2 Wl) . . . rg^w^ . . . rg m - l w n ) . (62) 

The code W describes the periodic orbit of the SB which is continued 
(f)(g) times to become periodic in the ST. Applying a generalized cyclic 
permutation on W is equivalent to applying the standard cyclic permutation 
on W. Applying a geometrical operation h on W is equivalent to operating 
letter by letter with h on W. The invariance group corresponding to W is 
H = C <E> Oh, where C is the group of cyclic permutations of order n ■ (f>(g). 
The simple decomposition of Ti is due to the commutativity of C and Oh, 
and it greatly facilitates the computations. 

3. If W is maximal with respect to the invariance group H, then the corre- 
sponding W is the representative of the periodic orbit. 

A comment on the uniqueness of this selection process is appropriate at this 
point. For any W we can uniquely construct the corresponding W and the 
invariance group and check the maximality of W. Hence, we are able to uniquely 
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decide whether W is a valid representative code or not. However, there are 
cases in which more than one W correspond to the same maximal W. It is 
straightforward to show that in these cases the basic code word W is symmetric 
under some operation(s): W = kW, k G Oh- To such symmetric codes must 
correspond symmetric periodic orbits, which is necessitated by the uniqueness 
theorem for the ST. But for the SB the symmetry of the orbit means that it is 
wholly contained in a symmetry plane, and therefore is not a proper classical 
orbit. Such orbits are nevertheless required for the semiclassical analysis and will 
be treated in the next section when dealing with semiclassical desymmetrization. 
In summary, we have shown so far that the mapping of a given proper periodic 
orbit to a code is well-defined and unique. 

In order for the coding to be useful and powerful, we need to establish unique- 
ness in the opposite direction, that is to show that for a given (unsymmetrical) 
W there corresponds at most one (proper) classical periodic orbit. The mapping 
W i— > W is very useful in that respect. Indeed, if there were two distinct periodic 
orbits of the SB with the same coding W, then we could repeat them <fi(g) times 
to get two distinct periodic orbits of the ST with the same code W, which is in 
contradiction with the theorem above. This proves the uniqueness of the relation 
between codes and periodic orbits. 

To facilitate the actual computation of periodic orbits of the SB, we have to 
establish their minimality property, similarly to the ST case. We need to prove 
that the length of a periodic orbit is a minimum, and that it is the only minimum. 
The minimality of a periodic orbit of the SB is proven by using again the unfolding 
to periodic orbits of ST, and noting that a minimum of is necessarily also a 
minimum of L^, since the latter is a constrained version of the former. Thus, 
periodic orbits of the SB are minima of L^. We finally have to show that there 
exists only a single minimum of L^. The complication here is that, in principle, 
a minimum of does not necessarily correspond to a minimum of L^, since 
there are, in general, more variables in the latter. We resolve this difficulty by 
using arguments from the proof of Schanz |HJ as follows. A necessary condition 
for minimality is that orbits are either externally reflected from the scatters or 
cut through them in straight segments. Internal reflections are not allowed for a 
minimum. Thus, if we extend a minimum of SB to ST, we necessarily get an orbit 



with no internal reflections. According to Schanz ||32|| , there is exactly one such 



orbit, which is the minimum in ST. This proves the uniqueness of the (global) 
minimum of in SB. 

These results allow us to use essentially the same algorithm as for the ST for 
the systematic search of periodic orbits of the SB. We need to extend the codes 
and the length functions to include a group element g, and to modify the rules 
according to which we choose an admissible and lexically maximal code word W. 
One also has to modify the computations of the monodromy matrix, as described 
in appendix [J]. 



49 



4.3 The properties and statistics of the set of periodic 
orbits 



The algorithm described above is capable of finding all of the periodic orbits up 
to any desired length. Before discussing the properties of this set, we find it 
appropriate to display a few typical periodic orbits, which were computed for the 
desymmetrized billiard with R = 0.2 (and S = 1.). The orbits are represented in 
an unfolded way in figures |16HT^- 

In this subsection we shall study in detail the spectrum of lengths of periodic 
orbits, a small interval thereof is shown in figure |2D|. Each horizontal strip provides 
the lower end of the length spectrum of Sinai billiards with 0.02 < R < 0.36. The 
spectrum corresponding to the lowest value of R shows clustering of the lengths 
near the typical distances of points of the Z 3 lattice (1, y/2, \/3, 2, . . .). Once 
R is increased, some of the periodic orbits which were allowed for the smaller 
R are decimated because of the increased effect of shadowing. However, their 
lengths become shorter, resulting in the proliferation of the periodic orbits with 
their length. This is best seen in the spectrum which corresponds to the largest 
value of R — the graphics is already not sufficiently fine to resolve the individual 
lengths. 

After these introductory comments, we now study the length spectrum in 
detail, and compare the theoretical expectations with the numerical results. The 
exponential proliferation of the periodic orbits puts a severe limit on the length 
range which we could access with our finite computer resources. However, we 
were able to compute the periodic orbits for a few values of the radius R, and 
concentrated on the R = 0.2 case in order to be able to perform a semiclassical 
analysis of the longest quantal spectrum (see next section). For this radius we 
found all the 586,965 periodic orbits up to length 5. This number of periodic 
orbits includes repetitions and time-reversed conjugates. We also computed for 
this radius all the 12,928,628 periodic orbits up to length 10 which have no more 
than 3 reflections. This comprises the database on which we based our further 
numerical studies and illustrations. The systematic algorithm which was used 
to produce this data set, together with a few tests which will be described here 
and in the next section, lead us to believe that the data set is both accurate and 
complete. 

Periodic orbits are expected to proliferate exponentially (e.g., 0). That is, 
the number N\ en (l) of periodic orbits of length less than I should approach asymp- 
totically §: 

AUO^^j^, l-oo, (63) 
where A is the topological entropy (per unit length). To examine the validity of 
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length = .600000 length = .770918 




length = 2.365230 length = 4.996244 




Figure 16: A sample of periodic orbits of the desymmetrized 3D Sinai billiard 
with S = 1, R = 0.2 with a single reflection. The periodic orbits are shown 
in the unfolded representation. The "full" spheres are those from which the 
periodic orbit reflects. The "faint" dotted spheres are those from which there is 
no reflection. 
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length = 4.236999 



Figure 17: A sample of periodic orbits of the 3D SB with 2 reflections. 
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length = 4.943742 




Figure 19: A sample of periodic orbits of the 3D SB with 7 reflections. The 
bottom periodic orbit undergoes 8 reflections. 
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Figure 20: Length spectra of periodic orbits for Sinai billiards with R values 
between 0.02 and 0.36 in steps of AR = 0.02. The vertical bars indicate the 
lengths of periodic orbits. 
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the above formula in our case we use the numerical data to compute: 




A„ um (/) = 7 ln V Lj , (64) 



where L eTg is a length below which we do not expect universality (i.e, the law 
to hold. The exponential proliferation implies: 

A num (/)^yln|e A '-e AL -|-^^A, I - oo . (65) 

Therefore, we expect A num (/) to approach a constant value A when I is sufficiently 
larger than L erg . In figure |2l] we show the results of the numerical computation 
of A num for the R = 0.2 database and for L erg = 2.5. The figure clearly indicates 
a good agreement between the data and the theory (|65D for A = 3.2. 

One of the hallmarks of classically ergodic systems is the balance between the 
proliferation of periodic orbits and their stability weights due to ergodic coverage 
of phase space. This is a manifestation of the uniform coverage of phase space 



and is frequently referred to as the "Hannay - Ozorio de Almeida sum rule" [[71 
It states that: 



PO 



where L p is the primitive length and Mj is the stability (monodromy) matrix 
(see appendix [J] for explicit expressions). The above relation is meaningful 
only after appropriate smoothing. For generic billiards the only classical length 
scale is the typical length traversed between reflections, and we expect (|66]) to 
approximately hold after a few reflections. In the Sinai billiard we are faced 
with the problem of an "infinite horizon", that is, that the length of free flight 
between consecutive reflections is unbounded. This is just another manifestation 



of the existence of the bouncing-ball families. According to |39], [5J this effect is 
responsible for a non-generic power-law tail in p(l): 

P(D « 1 - ^ , (67 ) 

where a(R) is a parameter that depends on the radius R. When R increases 
the influence (measure in configuration space) of the bouncing-balls is reduced, 
and we expect a(R) to decrease. To check ( |67|) we computed numerically the 
cumulant: 
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exact 
fit A=3.1 
fit A=3.2 
fit A=3.3 

2.5 3.0 3.5 4.0 4.5 5.0 

I 

Figure 21: The quantity A num (c.f. RHS of equation (0)) computed from the 
periodic orbit database of R = 0.2. We used L CTg = 2.5. The theoretical fit is 
according to equation (j65|) . 
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which should be compared to the theoretical expectation: 

P{1) = {l- L crg ) - a(R) In (-}-) . (69) 

The results are shown in figure 22. We considered R = 0.2 and 0.3 and included 
periodic orbits up to L max = 10 with number of reflections n < 3. The restriction 
on n facilitates the computation and is justified for moderate values of / since 
the contributions from higher n's are small. The observed deviation between 
the theoretical and numerical curves for R = 0.3 at / > 8 is due to the fact 
that periodic orbits with n = 4 become significant in this region. The above 
numerical tests confirm the validity of (|67|), with a(R) which is a decreasing 
function of R. In particular, for the length interval considered here, there is a 
significant deviation from the fully ergodic behavior (|66|). 

The sum-rule (|66|) which formed the basis of the previous analysis is an ex- 
pression of the ergodic nature of the billiards dynamics. In the next subsection we 
shall make use of similar sum-rules which manifest the ergodicity of the Poincare 
map obtained from the billiard flow by, e.g., taking the surface of the sphere and 
the tangent velocity vector as the Poincare section. The resulting return-map 
excludes the bouncing-ball manifolds since they do not intersect the section. 
However, their effect is noticed because between successive collisions with the 
sphere the trajectory may reflect off the planar faces of the billiard an arbitrary 
number of times. Thus, the number of periodic orbits which bounce n times off 
the sphere (n-periodic orbits of the map) is unlimited, and the topological en- 
tropy is not well defined. Moreover, the length spectrum of n-periodic orbits is 
not bounded. These peculiarities, together with the fact that the symbolic code 
of the map consists of an infinite number of symbols, are the manifestations of 
the infinite horizon of the unfolded Sinai billiard. The return map itself is dis- 
continuous but it remains area preserving, so the formulas which we use below, 
and which apply to generic maps, can be used here as well. 

The classical return probability is defined as the trace of the n-step classical 
evolution operator (see, e.g., |72] and references therein). It is given by: 

"<») - E iasfey, • (™> 

where n is the number of times the periodic orbit reflects from the sphere, V n is 
the set of all n-periodic orbits, n Pt j is the period of the primitive periodic orbit 
of which j is a repeated traversal. As a consequence of the ergodic nature of the 
map U(n) — > 1 in the limit n — > oo. However, due to the effect of the infinite 
horizon, the number of periodic orbits in V n is infinite, and in any numerical 
simulation it is important to check to what degree the available data set satisfies 
the sum rule. For this purpose we define the function: 

[iet( 7-M,)\ e «- L ^ < 71 > 

j inPn 
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R-0-3, n <3 

fully ergodic (slope 1) 

fitQ=2,l 
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Figure 22: The function P{1) (c.f. RHS of equation (j68|) ) computed for = 0.2 
and 0.3 and fitted according to equation d69|) . We also show the asymptotic 
prediction (IHO). 
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which takes into account only n-periodic orbits with Lj < I. In figure [23| we 
plot U(l; n) for R = 0.4 and n = 1, 2, 3. The results clearly indicate that for the 
present data saturation is reached, and once n > 2 the asymptotic value is very 
close to 1. Even at n — 1 one gets U(n = 1) w 0.8 which is surprisingly close 
to 1, bearing in mind that we are dealing with the fixed points of the map! It 
should be noted that to reach saturation in the case R = 0.4, n = 3 one needs 
536,379 periodic orbits up to I — 12, whose computation consumes already an 
appreciable amount of time. Thus, we are practically restricted to the few lowest 
n's in our computations. As can be seen in figure |23| the function dll(l; n)/dl is 
mostly supported on a finite interval of L values. Its width will be denoted by 
AL(n). 



4.4 Periodic orbit correlations 

In the previous subsection we discussed various aspects of the one-point statistics 
of the classical periodic orbits, and demonstrated their consistency with the stan- 
dard results of ergodic theory. Here, we shall probe the length spectrum further, 
and show that this spectrum is not Poissonian. Rather, there exist correlations 
between periodic orbits which have far-reaching effects on the semiclassical the- 
ory of spectral statistics of the quantum billiard. The semiclassical theory will 
be dealt with in section [7], and here we restrict ourselves to purely classical in- 
vestigations. 

Above we introduced the Poincare return map of the sphere, and have shown 
that the ergodicity of this map implies a sum rule for the set of n-periodic orbits 
of the map. We define the weighted density of lengths of n-periodic orbits as 
follows: 

d cl (l;n)= Y^A^l-Lj), (72) 



where Aj are given by: 



Aj \det(I-M 3 )\^ ' (73) 

where bj is the number of times the trajectory reflects from the planar boundaries. 
The amplitudes Aj are related to the standard semiclassical amplitudes Aj defined 
in (|9~9|) by Aj = nnajAj/Lj. 

The density ( [72] ) is different from the density p(l) defined previously (|BT)|) 
since: (a) it relates to the subset of the n-periodic orbits of the return map of 
the sphere, (b) it assigns a signed weight to each of the 5-functions located at a 
particular length, and (c) the absolute value of the weights in ( [72] ) are the square 
roots of the weights in (|66D. Densities with signed weights are not encountered 
frequently in spectral theory, but they emerge naturally in the present context. 
At this point the definition of d c \(l; n) might look unfamiliar and strange, but the 
reason for this particular choice will become clear in the sequel. 
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Figure 23: Upper plot: The function U(l;n) (c.f. equation (ffl|) ) for the cases 
i? = 0.4, n = 1,2,3. Lower plot: The function dU(l;n)/dl for the same cases. 
Both plots indicate the saturation of the classical return probability in spite of 
the infinitely many periodic orbits in V n . 
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To examine the possible existence of correlations in the length spectrum, we 
study the corresponding autocorrelation function: 



R d (Sl;n)= / dld d (l + 5l/2;n)d cl (l-5l/2;n) 



(74) 



The two-point form factor is the Fourier transform of R c i(5l;n), and it reads 
explicitly as: 

2 



K cl (k;n) 



+00 



e tkx R c i(x; n)dx 



Aj exp(ikLj] 



(75) 



The form factor has the following properties: 



K c i(k; n) is a Fourier transform of a distribution and therefore it displays 
fluctuations, which become stronger as the number of contributing orbits 
increases. Therefore, any discussion of this function requires some smooth- 
ing or averaging. We shall specify the smoothing we apply in the sequel. 

At k = 0, 

2 

K cl (0-n)= J2 A i ■ ( 76 ) 

Because of the large number of periodic orbits, the sum of the signed am- 
plitudes is effectively reduced due to mutual cancellations. Its value can be 
estimated by assuming that the signs are random. Hence, 

2 



£4 



(77) 



which will be shown below to be bounded. 
At large values of k, 

K cl (k;n) « ^ 9j A 



for k — > 00 



(7* 



where gj is the number of isometric periodic orbits of length Lj. Since large 
fluctuations are endemic to the form factor, this relation is meaningful when 
fc-averaging is applied. Comparing the last sum with (70) we can write: 

K c i(k; n) w (n p jgj) U(n) , for k — > 00 . (79) 

In our case of the 3D SB, n p j = n for the large majority of the periodic 
orbits in V n , which is the generic situation for chaotic systems. Also, gj = 2 
for almost all the periodic orbits with n > 3. Thus, one can safely replace 
( n p,j9j) whh In for large n. Moreover, as we saw above, U(n) —>■ 1 for large 
n, hence K c \ — > 2n for large k and n. 
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• If the length spectrum as defined above were constructed by a random 
sequence of lengths with the same smooth counting function U(l;n), or if 
the phases were picked at random, one would obtain the Poisson behavior 
of the form factor, namely, a constant: 

K d (k] n) « (n p g p )U(n) , for k > . (80) 

Here, AL(n) is the effective width of the length distribution defined above. 

Thus, we could identify two-point correlations in the classical length spectrum by 
computing K c \(k;n) and observing deviations from the fc-independent expression 



4.4.1 Numerical tests 

We used the periodic orbit database at our disposal to compute the form factors 
for several values of n and R. In each case presented we made sure that the 
function U (I; n) is numerically saturated. This guarantees that the (infinitely 
many) neglected periodic orbits have very small weight, and are thus insignificant. 
In figure [24] we present the numerical results, where we plotted the function: 

C cl {k; n) = 1 / dk' K cl (k'; n) , (81) 

™ ">min </fc min 

designed to smooth the fluctuations in K c \{k;n) |TT| . We started the integration 
at /c m j n > to avoid the large peak near k — 0, which otherwise overwhelms the 
results. In any case, the neglected small-fc region is irrelevant for the semiclassi- 
cal theory of quantal spectral correlations. Analyzing the results, we note that 
the asymptotic form factors (denoted as "Full classical form factor") approach 
constant values, which are indeed close to 2n, as predicted. More importantly, 
the deviations from the constant (Poissonian) result at low k demonstrate un- 
ambiguously the existence of correlations in the classical spectra. The structure 
of the form factor indicates that the classical spectrum is rigid on the scale of a 
correlation length A(n; R), which can be defined as the inverse of the k value at 



which the form factor makes its approach to the asymptotic value ||1 1|| . In the 
following we shall describe a few tests which prove that the observed correlations 
are real, and not a numerical artifact or a trivial consequence of the way in which 
the length spectral density is defined. 

The spectral density d c \(l;n) has an effective finite width AL(n) which was 
defined above. The fact that the lengths are constrained to this interval induces 
trivial correlations which appear on the scale AL(n), and we should check that 
this scale is sufficiently remote from the correlation scale \(n;R). To this end, 
and to show that the observed classical correlations are numerically significant, 
we scrambled the signs of the weights Aj by multiplying each of them with a 
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Figure 24: The averaged classical form factor C c \{k]n) (c.f. (|31~D) of the 3D SB 
for several values of n and i?. We also plot the averaged form factors with signs 
of the amplitudes scrambled, without cross-family terms, and with amplitudes 
averaged over family (length-correlation only). See text for details. 
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randomly chosen sign. We maintained, however, the time-reversal symmetry 
by multiplying conjugates by the same sign. The resulting form factors, shown 
in figure ^ (denoted as "Scrambled signs"), are consistent with the Poissonian 
value 2n for essentially all k values, and the difference between the scrambled 
and unscrambled data is large enough to add confidence to the existence of the 
classical correlations. This indicates also that the correlations are not due to the 
effective width of d c \(l;n), since both the scrambled and unscrambled data have 
the same effective width. 

On the other extreme, one might suspect that the classical correlations are due 
to rigidity on the scale of one mean spacing between lengths of periodic orbits. 
This is certainly not the case, since the typical mean length spacing for the cases 
shown in figure ^| is 10 -3 -lCT 4 , which implies a transition to the asymptotic 
value for much larger fc-values than observed. We therefore conclude, that the 
correlation length A(n; R) is much larger than the mean spacing between neigh- 
boring lengths. This is the reason why various studies of the length-spectrum 
statistics |)t], [73[] claimed that it is Poissonian. Indeed it is Poissonian on the 
scale of the mean spacing where these studies were conducted. The correlations 
become apparent on a very different (and much larger) scale, and there is no con- 
tradiction. The coexistence of a Poissonian behavior on the short length scales, 



and apparent rigidity on a larger scale was discussed and explained in |TT| . It was 
suggested there that a possible way to construct such a spectrum is to form it 
as a union of N 3> 1 statistically independent spectra, all having the same mean 
spacing A, and which show spectral rigidity on the scale of a single spacing. The 
combined spectrum with a mean spacing A/N will be Poissonian when tested on 
this scale, since the spectra are independent. However, the correlations on the 
scale A will persist in the combined spectrum. A simple example will illustrate 
this construction. Take a random (Poissonian) spectrum with a mean spacing 
1. Generate a shifted spectrum by adding A ^> 1 to each spectral point and 
combine the original and the shifted spectra to a single spectrum. On the scale 
1 the combined spectrum is Poissonian. However, the fact that each spectral 
point is (rigidly) accompanied by another one, a distance A apart, is a correlation 
which will be apparent at the scale A only. We use this picture in our attempt to 
propose a dynamical origin of the length correlations. 

4.4.2 The dynamical origin of the correlations 

As was already mentioned, the idea that periodic orbit correlations exist orig- 
inates from the quantum theory of spectral statistics which is based on trace 
formulas. The classical correlations are shown to be a manifestation of a fun- 



damental duality between the quantum and the classical descriptions ||, [Qfl . 
However, the effect is purely classical, and hence should be explained in classical 
terms, without any reference to the quantum mechanical analogue. The essential 
point is to find the classical origin of the partition of the periodic orbits to inde- 
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pendent and uncorrelated families, as was explained in the previous section. So 
far, all the attempts to find the classical roots of these correlations failed, and till 
now there is no universal theory which provides the classical foundations for the 
effect. For the Sinai billiard in 3D there seems to exist a physical-geometrical 
explanation, which is consistent with our data, and which is supported by further 
numerical tests. 

Consider the Sinai billiard with a sphere with a vanishingly small radius. In 
this case, all the periodic orbits which are encoded by words W built of the 
same letters Wi are isometric, independently of the ordering of the letters or 
the attached symmetry element g. This phenomenon can be clearly seen in the 
spectrum of lengths corresponding to R = 0.02 in figure In this case, it is 
clear that the spectrum of lengths is a union of "families" of periodic orbits, 
each family is characterized by a unique set of building blocks Wi, which are 
common to the family members. When the radius R increases and becomes 
comparable to the linear dimension of the billiard, the approximate isometry and 
the resulting correlations breaks down, and one should use a more refined and 
restrictive definition of a family. The aim is to find a partition to families which 
will restrict the membership in a family to the smallest set, without losing any 
of the correlation features. The most restrictive definition of a family in the 
present context will be to include all the periodic orbits which share the same 
W = (wi, W2, . . . , w n ) part of the code and have different admissible g symmetry 
elements. Words which are built of the same letters but in a different order define 
different families. Since there are 48 possible g's, each family consists of at most 
48 members and will be denoted by fl(W). It should also be noted that the signs 
of the weights Aj within a family do not change with R since they reflect the 
parity of g. The partition of the set of periodic orbits in families is not particular 
to just a few orbits, but rather, is valid for the entire set. This partition is the 
proposed source of the correlations that were observed in the form factor. This 
concept is illustrated in figure |25|, and graphic representations of two families are 
displayed in figure The most outstanding feature which emerges from figure 
(p6[) is that the orbits occupy a very narrow volume of phase-space throughout 
most of their length, and they fan out appreciably only at a single sphere. 

The above arguments suggest that the main source of correlations are the 
similarities of orbits within each family Q(W). To test this argument we per- 
formed a numerical experiment, in which we excluded the inter-family terms of 
the form factor, leaving only the intra-family terms. This excludes family-family 
correlations and maintains only correlations within the families. The results are 
shown in figure |24] (denoted by "Neglecting cross-family contributions"). The 
obvious observation is that the form factors were only slightly affected, proving 
that periodic-orbit correlations do not cross family lines! Thus, the main source 
of correlations is within the families Q(W). We mention that very similar results 
are obtained if inter-family sign randomization is applied instead of the exclusion 
of cross-terms. We note, that in most periodic orbit and its time-reversal 
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Figure 25: Two periodic orbits which are members of the same family of the 
quarter 2D SB. The two periodic orbits have the same W, hence they reflect off 
the same discs. But they correspond to two different symmetry elements, and 
hence are different. For simplicity the illustration is made for the 2D SB, but the 
same principle applies also to the 3D Sinai billiard. Left: Unfolded representation, 
right: Standard representation. 
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Figure 26: Two families of periodic orbits of the 3D SB, represented in 
unfolded representation of the SB. The faint spheres do not participate in 
code. 



68 



conjugate do not belong to the same family. Thus, neglecting the cross-family 
terms leads to partial breaking of time-reversal, which we compensated for by 
rectifying the intra-family form factor such that it will have the same asymptotic 
value as the full one. 

It is interesting to check whether the correlations are due to the lengths or due 
to the size of the amplitudes. To examine that, we not only neglected the cross- 
family terms, but also replaced the amplitudes A, within each family by constants 
multiplied by the original signs, such that the overall asymptotic contribution of 
the family does not change. The results are also plotted in figure |24| (denoted as 
"Length correlations only"). The resulting (rectified) form factors display slightly 
diminished correlations. However there is no doubt that almost all correlations 
still persist. This proves, that the correlations between the magnitudes of the 
weights play here a relatively minor role, and the correlations are primarily due 
to the lengths. 

There are a few points in order. First, the numerical results presented here 
concerning the classical correlations are similar to those of reference | JLljj . How- 
ever, here we considered the classical mapping rather than the flow, and this 
reduces the numerical fluctuations significantly. Using the mapping also enables 
the quantitative comparison to the semiclassical theory, which will be discussed 
in section [7|. Second, it is interesting to enquire whether the average number of 
family members iVf am (n;.R) increases or decreases with n. Since, if it decreases, 
our explanation of the origin of correlations becomes invalid for large n. The nu- 
merical results clearly indicate that iVf am (n; R), computed as a weighted average 
with the classical weights, increases with n, which is encouraging. For example, 
for the case R = 0.4 we obtained N iam = 9.64, 18.31, 21.09, 28.31 for n = 1, 2, 3, 
4, respectively. 

Thus, we were able to identify the grouping of orbits into "families" with the 
same code word W but with different symmetry g as the prominent source of the 
classical correlations in the 3D Sinai billiard. In each family the common geo- 
metric part of the code W sets the mean length and the different group elements 
g introduce the modulations. This pattern repeats for all the families, but the 
lengths of different families are not correlated. This finding conforms very well 
with the general scheme which was proposed to explain the typical correlations 
in the classical spectrum [f|, 74 1. However, a classical derivation of a quantitative 
expression for the correlations length A is yet to be done. 



4.4.3 Length correlations in the 3— Torus 

The ideas developed above about the correlations between periodic orbits in the 
fully chaotic billiard, have an analogue in the spectrum of lengths of periodic tori 
in the integrable case of the 3-torus. In section ^| we studied the quantum 3-torus 
of size S and showed in section that due to number theoretical degeneracies, 
the quantum form factor is not Poissonian. The form factor displays a negative 
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(repulsive) correlation which levels off at r* = I/7 = l/(2Sk). This can be 
transcribed into an expression for the correlation length of the classical spectrum 
in the following way. 

Expressing r* in units of length we obtain: 



2nd(k)r* = — . (82) 
2tt 



Consequently: 



from which we read off 



k*(L) = *£, (83) 

X(L) = ^ . (84) 

Since the lengths of the periodic orbits are of the form Lj = S x ^integer, the 
minimal spacing between periodic orbits near length L is: 

A min (L) = ^ , (85) 

and therefore 

A(L) = 2A min (L) . (86) 

In other words, the classical correlation length of the 3-torus coincides (up to 
a factor 2) with the minimal spacing between the periodic orbits. Therefore, 
A(L) indeed signifies the correlation length scale between periodic orbits, which 
is imposed by their number-theoretical structure. 
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5 Semiclassical analysis 



In the previous sections we accumulated information about the quantum spec- 
trum and about the periodic orbits of the 3D Sinai billiard. The stage is now 
set for a semiclassical analysis of the billiard. We shall focus on the analysis of 
the semiclassical Gutzwiller trace formula that reads in the case of the Sinai 
billiard: 

oo 

d(k) = ^25(k- k n ) « d{k) + d hh (k) + Aj cos(kLj) . (87) 

n=l PO 

The quantum spectral density on the LHS is expressed as the sum of three terms. 
The term d is the smooth density of states (see section |). The term G?t,b consists of 
the contributions of the non-generic bouncing-ball manifolds. It contains terms 
of the form (|3T| ) with different prefactors (which are possibly 0) due to partial or 
complete shadowing of the bouncing-ball family by the sphere. The last term is 
the contribution of the set of generic and unstable periodic orbits, where Lj denote 
their lengths and Aj are semiclassical amplitudes. One of the main objective 
of the present work was to study the accuracy of (|87|) by a direct numerical 
computation of the difference between its two sides. This cannot be done by a 
straightforward substitution, since three obstacles must be removed: 

• The spectrum of wavenumbers k n was computed for the fully desymmetrized 
Sinai billiard. To write the corresponding trace formula, we must remem- 
ber that the folding of the Sinai torus into the Sinai billiard introduces 
new types of periodic orbits due to the presence of symmetry planes, edges 
and corners. Strictly speaking, the classical dynamics of these orbits is 
singular, and becomes meaningful only if proper limits are taken. As ex- 
amples we mention periodic orbits that bounce off a corner, or that are 
wholly confined to the symmetry planes. These periodic orbits are isolated 
and unstable, and should not be confused with the bouncing-ball families 
which are present both in the ST and in the SB. For periodic orbits that 
reflect from a corner but are not confined to symmetry planes, the diffi- 
culty is resolved by unfolding the dynamics from the SB to the ST as was 
described in the previous section. Periodic orbits which are confined to 
symmetry planes are more troublesome since there is more than one code 
word W which correspond to the same periodic orbit. We denote the latter 
as "improper". The 3D Sinai billiard is abundant with improper periodic 
orbits, and we cannot afford treating them individually as was done e.g. by 
Sieber |5D| for the 2D hyperbola billiard. Rather, we have to find a general 



and systematic method to identify them and to calculate their semiclassical 
contributions. This will be done in the next subsection. (The semiclassi- 
cal contributions of the improper periodic manifolds for the integrable case 
R = were discussed in section 
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As it stands, equation fl8~7D is a relation between distributions rather than 
between functions, and hence must be regulated when dealing with actual 
computations. Moreover, even though our quantum and classical databases 
are rather extensive, the sums on the two sides of the equation can never be 
exhausted. We overcome these problems by studying the weighted "length 
spectrum" obtained from the trace formula by a proper smoothing and 
Fourier-transformation. It is defined in subsection B.2L 



Finally, we must find ways to rid ourselves from the large, yet non-generic 
contributions of the bouncing-ball families. This was achieved using rather 



elegant tricks which are described in subsections |5.4|, |5.5| below 



5.1 Semiclassical desymmetrization 



To derive the spectral density of the desymmetrized Sinai billiard we make use 
of its expression in terms of the (imaginary part of the) trace of the SB Green 
function. This Green function satisfies the prescribed boundary conditions on all 
the boundaries of the fundamental domain, and the trace is taken over its volume. 
In the following we shall show how to transform this object into a trace over the 
volume of the entire ST, for which all periodic orbits are proper (no symmetry 
planes). This will eliminate the difficulty of treating the improper orbits. To 



achieve this goal we shall use group-theoretical arguments [[H], |6(], [75], [76], [77]. 



The final result is essentially contained in [78 



When desymmetrizing the ST into SB, we have to choose one of the irreps of 
Oh to which the eigenfunctions of the SB belong (see section |2.2.3|). We denote 



this irrep by 7. We are interested in the trace of the Green function of the SB 
over the volume of the SB which is essentially the density of states: 



T = Tr SR G 



(7), 



SB^SB 



r, r 



(88) 



One can apply the projection operation 
function of the ST: 



and express Ggg using the Green 



G 



(7), 

SB 



r, r 



7 9&O h 



(89) 



where X {$) is the character of g in the irrep 7 and Z 7 is the dimension of 
7. It can be verified that the above Gsb satisfies the inhomogeneous Helmholtz 
equation with the correct normalization, and it is composed only of eigenfunction 
that transform according to 7. Thus: 



T =T Ex (7) *^)Tr S BG 

7 g&o h 



ST[r,gr 



(90) 
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To relate TrsB with Tr$T we use the relation: 



G ST (f,f') = G ST (y,hf') VheO h (91) 

which can be proven by e.g. using the spectral representation of Gst- In partic- 
ular, we can write: 

heo h 

Combining fl9"2] ) with ( pOj ) we get: 

T = i^T E X {7) *(^)TrsBGsT(^r,^f') 

= 4^T E X (7) *(^ 1 )Trs B Gs T (/ i r ) (^/ i - 1 )/ i f') 

= i^T E X (7) *(^)TrsBG S T(/ i r,^f / ). (93) 

7 h.fceOh 

To obtain the second line from the first one, we recall that the character is the 
trace of the irrep matrix, and we have in general Tr(ABC) = Tr(CAB), therefore 
x{g) — x{hgh~ l ). The third line is obtained from the second one by fixing h and 
summing over g. Since hgih^ 1 = h^h^ 1 g\ = §2 the summation over g is a 
rearrangement of the group. We now apply the geometrical identity: 

[ d 3 rf(hr) = I d 3 r/(f) (94) 

hdO h 



SB ./ST 



to cast (|93|) into the desired form: 

T = E X (7) *(^)Trs T GsT(r,^') , (95) 
7 geO h 

where we relabelled k as g for convenience. The result (|95|) is the desired one, 
since T is now expressed using traces over ST which involve no symmetry planes. 
Semiclassically, the formula (|95|) means that we should consider all the periodic 
orbits of the ST modulo a symmetry element g to get the density of states of the 
SB. Therefore, the difficulty of handling improper orbits is eliminated, since in 
the ST all of the isolated periodic orbits are proper. 

Let us elaborate further on ( P5| ) and consider the various contributions to it. 
A proper periodic orbit of the SB with code (W; g) has 48 realizations in the ST 
which are geometrically distinct. They are obtained from each other by applying 
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the operations of Oh- These conjugate periodic orbits are all related to the same 
g and thus have the same lengths and monodromies. Consequently they all have 
the same semiclassical contributions. Hence, their semiclassical contribution to T 
is the same as we would get from naively applying the Gutzwiller trace formula to 
the SB, considering only proper periodic orbits. This result is consistent with our 
findings about classical desymmetrization (section |4]^ above). For the improper 
periodic orbits there is a difference, however. There are genuine semiclassical 
effects due to desymmetrization for unstable periodic orbits that are confined 
to planes or to edges, notably large reduction in the contributions for Dirichlet 
conditions on the symmetry planes. 

To demonstrate this point, let us consider in some detail an example of the 
periodic orbit that traverses along the 8-fold edge AE in figure [I]. For the ST (no 
desymmetrization) its semiclassical contribution is: 

A, = £ . (96) 

271 

For the SB there are 8 code words that correspond to the periodic orbit (s) which 
traverses along this 8-fold edge. A calculation yields for the semiclassical contri- 
bution: 

R 



A 8 = 1 - 



2 ± 2Vi=2jj±/>(^f) 



(97) 



where (3 = R/S. The upper sign is for the case of the totally symmetric irrep, 
and the lower one for the totally antisymmetric irrep. In the antisymmetric case 
we get for (3 -C 1: 

which means that the desymmetrization greatly reduces the contribution of this 
periodic orbit in case of Dirichlet boundary conditions on the planes. For the case 
of our longest spectrum (R = 0.2, S — 1) this reduction factor is approximately 
2 x 10 -4 which makes the detection of this periodic orbit practically impossible. 
For Neumann boundary conditions the contribution is comparable to the ST case 
and is appreciable. 

The formula (f55|) together with the algorithm described above are the basis 
for our computations of the semiclassical contributions of the periodic orbits of 
the SB. Specifically, the contribution of a code W is given by: 



ir^r |det (I - M w 



where is the length of the periodic orbit, = (# of distinct realizations of 

W under O/J/48 and r is the repetition index. The term is due to the reflec- 
tions from the spheres and is determined by the boundary conditions on them. 
For Neumann boundary conditions = 1, for Dirichlet boundary conditions 
aw = (—1)™, where n is the number of bounces. 
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5.2 Length spectrum 

Having derived the explicit expression for the semiclassical amplitudes for the SB 
fl99|), we are in position to transform the trace formula fl87D to a form which can 
be used for numerical computations which test its validity. We define the length 
spectrum as the Fourier transform of the density of states: 

1 f +QC 1 
D(l) = —= I d(k) e ikl dk = —= J2 e lknl • (100) 

V 27T J_oo V 27T 



For convenience we define d{— k) = d(k) k_ n = —k n and the sum is carried 
out for all n G Z\{0}. Using the trace formula fl5?| ) we obtain semiclassically: 

D sc (l) = D(l) + D hh (l) + \J\Aj W ~ L 3 ) + 5(1 + Lj)] . (101) 

po ' 

In the above D(l) is a singularity at I = which is due to the smooth density 
of states. The length spectrum is sharply peaked near lengths of periodic orbits 
hence its name. To regularize ( |101| ) such that it can be applied to finite samples 
of the quantum spectrum, we use a weight function and construct the weighted 
length spectrum []73[| : 

1 r+°° i 
D^ w \l-k) = —= w(k-k')d(k')e ik ' l dk' = —=^w{k-k n ) e lknl (102) 

V 27T J -co V 27T 



where w is a weight function (with an effective finite support) that is concentrated 
at the origin. The corresponding semiclassical expression is: 

Di:\i ; k) = m(i) + d^\i) + H l - L^ k{l ~ Li) + *>{l + L 3 y k{ - l+L ^] , 

po 2 

(103) 

where w(l) = (l/y/2n) f_°° w(k) e lkl dk is the Fourier transform of w(k). 

In principle, d(k) and D(l) contain the same information and are therefore 
equivalent. However, for our purposes, it is advantageous to use the length spec- 
trum D(l) (and in practice D^"\l\ k)) rather than the spectral density fl87f) for 
the following reasons: 

• The regularized semiclassical length spectrum, is absolutely conver- 

gent for suitably chosen windows (e.g. Gaussians). This is in contrast 
with the original trace formula (|37j). 



There is an exact mathematical result |80[ that states that for billiards the 
singular supports of D(l) and of D sc (l) are the same, if the infinite spectra 
are considered. This exact quantum-classical result specifically relates to 
the length spectra. It is therefore useful to identify and treat transient 
effects (e.g. diffraction contributions) for finite spectra using the length 
coordinate. 
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The trace formula fl8T|) can be considered as a means to quantize a chaotic 
system, since it expresses the quantal density of states in terms of the clas- 
sical length spectrum. However, in practice this is not convenient because 
the semiclassical amplitudes are only leading terms in asymptotic series in 
k (equivalently in h). For finite values of k there can be large deviations 
due to sub-leading corrections || |7|] and also due to significant diffraction 



corrections || |HJ, p8|| . Treating the trace formula the other way ("inverse 
quantum chaology") is advantageous because the quantal amplitudes have 
all equal weights 1. 

The appearance of peaks in both d(k) and D(l) comes as a result of the 
constructive interference of many oscillatory contributions. Any missing 
or spurious contribution can blur the peaks (see figure |27] for an example 



with a single energy level missing). For the energy levels we have a good 
control on the completeness of the spectrum due to Weyl's law (see section 
|). As discussed above, this is not the case for periodic orbits where we 
do not have an independent verification of their completeness. Hence it is 
advantageous to use the energy levels which are known to be complete in 
order to reproduce peaks that correspond to the periodic orbits. 

For the Sinai billiard the low-lying domain of the spectrum is peculiar due 



to effects of desymmetrization (see section |2.4j) . For Dirichlet boundary 
conditions on the planes, the levels k n R < 9 are very similar to those of the 
integrable case (R = 0). The "chaotic" levels for which the semiclassical 
approximation is valid (k n R > 9) thus start higher up, which makes the 
semiclassical reproduction of them very difficult in practice even with the 
use of Berry-Keating resummation techniques ||69|| . On the other hand, 



using the quantum levels we can reproduce a few isolated length peaks, as 
will be seen in the sequel. 

In the following we shall demonstrate a stringent test of the completeness and 
of the accuracy of the quantal spectrum using the length spectrum. Then we 
shall investigate the agreement between the quantal and the semiclassical length 
spectra. We shall employ a technique to filter the effects of the bouncing balls, 
such that only generic contributions remain. 

5.3 A semiclassical test of the quantal spectrum 

In the following we use the length spectrum in order to develop a stringent test 
of the completeness and integrity of the quantal spectrum. This supplements the 
integrity and completeness analysis of the quantal spectrum done in subsections 
2~3)- |2.5| . The idea is to focus on an isolated contribution to the length spectrum 



that can be compared to an analytical result. In section 3.1 we discussed the 



integrable billiard (R = 0) and observed that there are contributions to the 
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density of states due to isolated but neutral periodic orbits. The shortest periodic 
orbit of this kind has length S/y/3 ~ 0.577S" and was shown in figure |6|. Its 
contribution must prevail for R > until it is shadowed by the inscribed sphere, 
which occurs at R = S/ ~ 0.41S*. Being the shortest bouncing ball, its length 
is distant from the other bouncing balls. As for the isolation from other generic 
periodic orbits, for R = 0.2 there is a nearby contribution of the shortest unstable 
periodic orbit of length 0.65*. However, for Dirichlet boundary conditions on the 
planes the latter is practically eliminated due to symmetry effects as was discussed 



in section |5.1| . Since other periodic orbits are fairly distant, this shortest bouncing 
ball is an ideal test-ground of the length spectrum. Using and a Gaussian 
window: 



w(k - k') 



1 



V2 



exp 



(k - k'f 
2a 2 



(104) 



one obtains the contribution of the shortest bouncing ball to the length spectrum: 



D 



(«0 

sc, shortest— bb 



(l;k) 



Jk(l-S/VS) 

(6tt) 3 /2 



exp 



-{I- S/Vs) 2 a 2 /2 



(105) 



Due to its isolation, one expects that the shortest bouncing ball gives the dom- 
inant contribution to the length spectrum near its length. Thus, for I « S/ y/3, 



one has 



I n( w )| 



\D 



sc, shortest— bb I 



The latter is independent of A;. To test the 



above relation, we computed the quantal length spectrum for R = and 
R = 0.2 for two different values of k, and compared with (|105|) . The results are 
shown in figure |27|, and the agreement is very satisfactory. 

To show how sensitive and stringent this test is, we removed from the R = 0.2 
quantal spectrum a single level, ki 500 = 175.1182, and studied the effect on the 
length spectrum. As is seen in the figure, this is enough to severely damage the 
agreement between the quantum data and the theoretical expectation. Therefore 
we conclude that our spectrum is complete and also accurate to a high degree. 



5.4 Filtering the bouncing-balls I: Dirichlet— Neumann dif- 
ference 

The final goal of our semiclassical analysis is to test the predictions due to 
Gutzwiller's trace formula. Since the 3D Sinai is meant to be a paradigm for 
3D systems, we must remove the influence of the non-generic bouncing-ball fam- 
ilies and find a way to focus on the contributions of the generic and unstable 
periodic orbit. This is imperative, because in the 3D Sinai billiard the bouncing 
balls have contributions which are much larger than those of the generic peri- 
odic orbits. Inspecting equations ( j3lf) and (0), we find that the contributions of 
the leading-order bouncing balls are stronger by a factor of k than those of the 
generic periodic orbits. This is worse than in the 2D case, where the factor is \fk. 
To show how overwhelming is the effect of the bouncing balls, we plot in figure |2S] 
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Figure 27: Absolute value of the quantal length spectra |Z)W| with a Gaussian 
window, a = 30, compared to the theoretical prediction (|105|) . The location of 
the shortest bouncing ball is indicated by the vertical line. 
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the quantal lengths spectra for R = and R = 0.2 (Dirichlet everywhere) 

together with |-D^ | which contains contributions only from generic and unstable 
periodic orbits. One observes that all the peaks in the quantal length spectra 
are near lengths of the bouncing balls. Contributions of generic periodic orbits 
are completely overwhelmed by those of the bouncing balls and cannot be traced 
in the quantal length spectrum of R = 0.2. Also, we see that for R = 0.2 the 
peaks are in general lower than for R = 0. This is because of the (partial or 
complete) shadowing effect of the inscribed sphere that reduces the prefactors of 
the bouncing balls as R increases. 

In the case of the 2D Sinai billiard it was possible to analytically filter the 



effect of the bouncing balls from the semiclassical density of states |37| , I7| . In 
three dimensions this is much more difficult. The functional forms of the contri- 
butions to the density of states of the bouncing balls are given in (pT|) , but it is 
a difficult geometric problem to calculate the prefactors which are proportional 
to the cross sections of the bouncing-ball manifolds in configuration space. The 
desymmetrization makes this difficulty even greater and the calculations become 
very intricate. In addition, there is always an infinite number of bouncing-ball 
manifolds in the 3D Sinai. This is in contrast with the 2D Sinai, in which a finite 
(and usually quite small for moderate radii) number of bouncing-ball families 
exist. All this means, that an analytical subtraction of the bouncing-ball contri- 
butions is very intricate and vulnerable to errors which are difficult to detect and 
can have a devastating effect on the quantal-semiclassical agreement. 

In order to circumvent these difficulties, we present in the following an efficient 
and simple method to get rid of the bouncing balls. The idea is simple: The 
bouncing balls are exactly those periodic orbits that do not reflect from the 
sphere. Therefore, changing the boundary conditions on the sphere does not 
affect the bouncing-ball contributions. Thus, the semiclassical density of states 
for Dirichlet / Neumann boundary conditions on the sphere is: 

d D / N = Jd/n + d hh + dj^f . (106) 

The difference dn — d^ is hence independent (in leading approximation in k) of 
G?bb and has the standard form of a trace formula: 

d D _ N = d D {k) - d^{k) = [d B {k) - d N {k)] + ( A ? } ~ A T ] ) cos(ArLj) . (107) 

PO 

Here Af\ Af ] are the coefficients that correspond to Dirichlet and Neumann 
cases, respectively. In fact, for Dirichlet, each reflection with the sphere causes a 
sign change, while for Neumann there are no sign changes. Therefore: 

■ (d-n) _ a(d) _ ,(N) _ f 2Aj odd number of reflections (108) 
3 3 3 I even number of reflections 
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Figure 28: Quantal length spectra for R = and R = 0.2 compared to semiclas- 
sical length spectrum for R = 0.2 that contains only generic, unstable periodic 
orbits. In all cases k = 160, a = 30. The locations of the bouncing balls are 
indicated: Daggers for 2-parameters bouncing balls that occupy 3D volume in 
configuration space, stars for 2D bouncing balls and crosses for ID bouncing 
balls. 
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and we expect to observe in the length spectrum of g?d-n contributions only due 
to generic periodic orbits with an odd number of reflections. The results of the 
numerical computations are presented in figure ^ where we compare the quantal 
(exact) vs. semiclassical (theoretical) length spectra. We observe on the outset 
that in contrast to figure the quantal and semiclassical length spectra are of 
similar magnitudes and the bouncing balls no longer dominate. The peaks near 
lengths that correspond to the bouncing balls are greatly diminished, and in fact 
we see that the peak corresponding to the shortest bouncing ball (I ~ 0.577) 
is completely absent, as predicted by the theory. Even more important is the 
remarkable agreement between the quantal and the semiclassical length spectra 
which one observes near various peaks (e.g. near I = 0.75, 1.25, 2). Since the 
semiclassical length spectrum contains only generic contributions from unstable 
periodic orbits, this means that we demonstrated the existence and the correct- 
ness of these Gutzwiller contributions in the quantal levels. Therefore, one can say 
that at least as far as length spectra are concerned, the semiclassical trace formula 
is partially successful. There are, however, a few locations for which there is no 
agreement between the quantal and the semiclassical length spectra. The places 
where this discrepancy takes place are notably located near 3D bouncing-ball 
lengths. This suggest that there are "remnants" of the bouncing-ball contribu- 
tions that are not filtered by the Dirichlet - Neumann difference procedure. It 
is natural to expect that these remnants are most prominent for the strongest 
(3D) bouncing balls. The origin of these remnants are the periodic orbits that 
are exactly tangent to the sphere. As an example, consider the 3D bouncing-ball 
families that are shown in figure ^ (upper part). The tangent orbits that are 
related to them constitute a 1-parameter family that surrounds the sphere like 
a "corona". For a single tangent traversal their contributions acquire opposite 
signs for Dirichlet and Neumann boundary conditions on the sphere. Hence the 
Dirichlet - Neumann difference procedure still include these contributions which 
is apparent in the large discrepancy near 1 = 1. For two tangent traversals the 
Dirichlet and Neumann contributions have the same sign and hence cancel each 
other. This is indeed confirmed in figure |29| where we observe that near / = 2 
there is no discrepancy between the quantal and the semiclassical length spectra. 

The above mentioned tangent orbits belong to the set of points in phase 
space in which the classical mapping is discontinuous. Semiclassically they give 
rise to diffraction effects. Tangent orbits were treated for the 2D case in our work 



H, |3q| . To eliminate their effects we hence need to sharpen our tools and to 
find a better filtering method than the present Dirichlet - Neumann difference 
procedure. This is performed in the following using mixed boundary conditions. 
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Figure 29: Dirichlet-Neumann difference length spectra for R = 0.2, with k = 
100, a = 30. The semiclassical length spectrum is computed according to ( |108| ). 
The daggers, stars and crosses indicate the positions of the bouncing balls (refer 



to figure |28|) and the vertical bars indicate the positions of the generic, unstable 
periodic orbits. 
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5.5 Filtering the bouncing-balls II: Mixed boundary con- 
ditions 



The idea behind the Dirichlet-Neumann difference method was to subtract two 
spectra which differ only by their boundary conditions on the sphere. This can 
be generalized, if one replaces the discrete "parameter" of Dirichlet or Neumann 
conditions by a continuous parameter a, and studies the differences of the corre- 



sponding densities of states d(k; a\) — d(k; 0:2) ■ In section [2.1| we discussed the 
mixed boundary conditions regarding the exact quantization of the 3D SB and 
gave the a-dependent expressions for the quantal phase shifts. Mixed boundary 
conditions were extensively discussed in [^2], j33f . 

To include the mixed boundary conditions in the semiclassical trace formula 
we generalize the results of Berry ||15|| . There, he derived the trace formula for the 
2D Sinai billiard from an expansion of the KKR determinant in terms of traces. 
If one uses the 3D KKR matrix with (Q) and perform a similar expansion, the 
result is a modification of the Gutzwiller terms as follows: 

Aj cos(kLj) — > Aj cos(kLj + rijn + <pj) , (109) 



, k cot a . 
arctan | ^ ] . (110J 



»=i 



k cos 9 - 



Here Aj are the semiclassical coefficients for the Dirichlet conditions on the sphere 
(c.f. equation (|9~9])) and rij counts the number of reflections from the sphere. The 
angles #| are the reflection angles from the sphere measured from the normal of 
the j'th periodic orbit. It is instructive to note that the phases (fTTH) above are 
exactly the same as those obtained by a plane wave that reflects from an infinite 
plane with mixed boundary conditions @. This is consistent with the local 
nature of the semiclassical approximation. A prominent feature of the mixed 
boundary conditions which is manifest in (|109|) is that they do not affect the 
geometrical properties (length, stability) of the periodic orbits. Rather, they only 
cause a change of a phase which depends on the geometry of the periodic orbit. 
This is due to the fact that the mixing parameter a has no classical analogue. 
The invariance of periodic orbits with respect to a renders the mixed boundary 
conditions an attractive parameter for e.g. investigations of parametric statistics. 
This was discussed and demonstrated in detail in p3| . 



We are now in a position to apply the mixed boundary conditions to get an 
efficient filtering of the bouncing-ball contributions. We first note that if we fix 
K, then the levels are functions of a: k n = k n (a). Let us consider the derivative 
of the quantal counting function at a = 0: 



da 



£ [A - A*(a)] 



da 

a=0 



a=0 



S3 



E 



dk n (a) 
do; 



5(k-k n ), (111) 

a=0 



where k n = k n (0) are the Dirichlet eigenvalues. Hence, the quantity d is a 
weighted density of states with delta-peaks located on the Dirichlet eigenvalues. 

The semiclassical expression for d does not contain the leading contribution of 
the bouncing balls, since this contribution is independent of a. The semiclassical 
contributions of the isolated periodic orbits to d are of the form AjBj cos(kLj), 
where 

5 i = ^5>s#. (112) 

i=l 

This is easily derived from ( |109| ) and ( |110| ). Since the reflection angles 9^ are 
in the range [0, 7r/2], the coefficient Bj vanish if and only if 9\ = rr/2 for all 
% = 1, . . . , rij, which is an exact tangency. Therefore, exactly tangent periodic 
orbits are also eliminated by the derivative method. This is the desired effect of 
the mixed boundary conditions method that serves to further clean the spectrum 
from sub-leading contributions of the bouncing balls. We summarize equations 
( |TTlD and (|IT2D : 



d(k) = V ^( k ~ K) W ( S ™°^ h ) + Yl A 3 B i CO < kL i) » V n = [ ) 

n ^ ^ PO ^ ' a=0 

(113) 

To check the utility of d and to verify ( |113|) we computed both sides of (|1 13|) for 
R = 0.2 and k = 100. The quantal spectrum was computed for a = 0.003 and the 
derivatives v n were obtained by the finite differences from the a = (Dirichlet) 
spectrum. The coefficients Bj were extracted from the geometry of the periodic 



orbits. In figures |30] and [Jl] the length spectra are compared. The agreement 
between the quantal and the semiclassical data is impressive, especially for the 
lower /-values. There are no significant remnants of peaks near the bouncing-ball 
locations, and the peaks correspond to the generic and unstable periodic orbits. 
This demonstrates the utility of using d as an efficient means for filtering the 
spectrum from the non-generic effects. 

The quantal-semiclassical agreement of the length spectra is not perfect, how- 
ever, and it is instructive to list possible causes of this disagreement. We first 
recall that the semiclassical amplitudes A, are the leading terms in an asymptotic 
series, hence we expect corrections of order 1/k to the weights of periodic orbits. 
They are denoted as h corrections and were treated in detail by Gaspard and 
Alonso and by Alonso and Gaspard 0. In our case, however, 1/k w 1/100 and 
these corrections are not expected to be dominant. More important are diffrac- 
tion corrections which are also finite k effects that stem from the existence of a 
concave component (the sphere) in the billiard. Several kinds of diffraction cor- 
rections to the trace formula were analyzed for 2D billiards. Vattay, Wirzba and 
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Figure 30: Length spectra for the mixed boundary conditions derivative method 
( |113| ). Data are for R = 0.2, k = 150, a = 30, k = 100. The dashed line represents 
\D^ w \l) — Dsc\l)\. The daggers, stars and crosses indicate the positions of the 
bouncing balls (refer to figure ^8|) and the vertical bars indicate the positions of 
the generic, unstable periodic orbits. 
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Figure 31: Continuation of figure [30] to 2.5 < / < 5. We did not indicate the 
locations of unstable periodic orbits due to their enormous density. 



Rosenqvist || considered creeping orbits, and we considered in |38| penum- 
bra corrections. (The penumbra is the region in phase space which is close to 
tangency: \£ — kR\ ~ (kR) 1 ^ 3 , where hi is the angular momentum.) We list the 
various diffraction corrections in the following: 

Creeping orbits: These are orbits which are classically forbidden. They "creep" 
over concave parts of the billiard, and their semiclassical contribution is ex- 
ponentially small in A; 1 / 3 , which should be negligible for the k values that 
we consider. 

Exactly tangent orbits: These were already mentioned above, and we showed 
that their contributions are eliminated to a large extent by the mixed 
boundary conditions procedure. For 2D systems we found, however, that 
this is true in leading order only, and there are small remnants of the tan- 
gent orbits in the weighted density d [[38] . The magnitude of the remnants 
in 2D is 0(l/yk), which is smaller than O(k ) of a generic unstable peri- 
odic orbit. In 3D, a similar analysis shows that the remnants of each family 
of tangent orbits is 0{k ) which is the same magnitude as for unstable 



periodic orbits. Reviewing figures 30 and 31, we can observe some of the 



peaks of the quantal-semiclassical difference near lengths that correspond 
to exactly tangent orbits. 

Unstable and isolated periodic orbits that traverse the penumbra: We 



have shown in [8~1, 133 that for periodic orbits which just miss tangency with 
a concave component of the billiard boundary, there is a correction to the 
semiclassical amplitude A, which is of the same magnitude as Aj itself. 
These 0(1) diffraction corrections are the most important corrections to 
the trace formula for generic billiards. For periodic orbits which reflect at 
an extreme forward direction from a concave component, the amplitude Aj 
is very small due to the extreme classical instability. If we include diffraction 
corrections, the semiclassical contributions of these orbits get much larger. 
Therefore, the semiclassical contributions of periodic orbit that traverse 
the penumbra must be radically corrected. Moreover, we found that if one 
considers all the periodic orbits up to the Heisenberg length Lh = 2ird(k) 
(necessary to obtain a resolution of one mean level spacing), then almost all 
of the periodic orbits are vulnerable to penumbra diffraction corrections. 

Classically forbidden periodic orbits that traverse the shaded part of 
the penumbra: Penumbra diffraction effects lead to semiclassical contri- 
butions from periodic orbits that slightly traverse through a concave com- 
ponent. Since they do not relate to classically allowed orbits, they represent 
new contributions to the trace formula rather than corrections of existing 
ones. Their magnitudes are comparable to those of generic unstable peri- 
odic orbits. 
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The above list of corrections, which was compiled according to studies of 2D 
billiards, suggests that there is a wealth of effects that must be considered if one 
wishes to go beyond the Gutzwiller trace formula. It is very difficult to implement 
these corrections systematically even for 2D billiards, and it goes beyond the scope 
of the present work to study them further for the 3D Sinai billiard. We mention 
in passing that except exact tangency, the penumbra effects are transient and 
depend on k. 

According to the mathematical theorem mentioned above, the quantal 
and the semiclassical length spectra are asymptotically the same. The signifi- 
cance of our findings in this section is that we have shown that the quantal- 
semiclassical agreement is achieved already for finite and moderate values of k 
and that the corrections are not very large (for the /-range we looked at). This 
is very encouraging, and justifies an optimistic attitude to the validity of the 
semiclassical approximation in 3D systems. However, obtaining accurate energy 
levels from the trace formula involves many contributions from a large number 
of periodic orbits. Therefore, one cannot directly infer at this stage from the 
accuracy of the peaks of the length spectrum to the accuracy of energy levels. 
There is a need to quantify the semiclassical error and to express it in a way 
which makes use of the above semiclassical analysis. This is done in section 
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6 The accuracy of the semiclassical energy spec- 
trum 



One of the most important applications of the trace formula is to explain the 
spectral statistics and their relation to the universal predictions of Random Ma- 
trix Theory (RMT) j|, |J. However, a prerequisite for the use of the semiclassical 
approximation to compute short-range statistics is that it is able to reproduce the 
exact spectrum within an error comparable to or less than the mean level spacing! 
This is a demanding requirement, and quite often it is doubted that the semi- 
classical approximation is able to reproduce precise levels for high-dimensional 
systems on the following grounds. The mean level spacing depends on the dimen- 
sionality (number of freedoms) of the system, and it is 0{h d ) Gutzwiller 
quotes an argument by Pauli |[7]] to show that in general the error margin for the 
semiclassical approximation scales as 0(h 2 ) independently of the dimensionality. 
Applied to the trace formula, the expected error in units of the mean spacing, 
which is the figure of merit in the present context, is therefore expected to be 
0(h 2 ~ d ). We shall refer to this as the "traditional estimate". It sets d = 2 as 
a critical dimension for the applicability of the semiclassical trace formula and 
hence for the validity of the conclusions which are drawn from it. The few sys- 
tems in d > 2 dimensions which were numerically investigated display spectral 
statistics which adhere to the predictions of RMT as accurately as their counter- 



parts in d = 2 ||53| , [54] , |19|| . Thus, the traditional estimate cannot be correct in 
the present context, and we shall explain the reasons why it is inadequate. 

In this section we shall develop measures for the accuracy of the semiclassical 
energy levels. We shall then derive formulas to evaluate these measures. Using our 
quantal and classical (periodic orbits) databases for the 2D and 3D Sinai billiards, 
we shall apply the formulas and get numerical bounds for the semiclassical errors. 

The problem of the accuracy of the energy spectrum derived from the semi- 
classical trace formula was hardly discussed in the literature. Gutzwiller quotes 
the traditional estimate of 0{h 2 ~ d ) [Q, 173]. Gaspard and Alonso M, Alonso and 
Gaspard @ and Vattay, Wirzba and Rosenqvist || derived explicit and generic 
h corrections for the trace formula, but do not address directly the issue of semi- 
classical accuracy of energy levels. Boasman [88] estimates the accuracy of the 
Boundary Integral Method (BIM) [|14]] for 2D billiards in the case that the exact 
kernel is replaced by its asymptotic approximation. He finds that the result- 
ing error is of the same magnitude as the mean spacing, in agreement with the 
traditional estimate. However, the dependence of the semiclassical error on the 
dimensionality is not established. We also mention a recent work by Dahlqvist 
J55fl in which the semiclassical error due to penumbra (diffraction) effects is ana- 
lytically estimated for the 2D Sinai billiard. The results are compatible with the 
ones reported here. 
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6.1 Measures of the semiclassical error 

In order to define a proper error measure for the semiclassical approximation of 
the energy spectrum one has to clarify a few issues. In contrast with the EBK 
quantization which gives an explicit formula for the spectrum, the semiclassical 
spectrum for chaotic systems is implicit in the trace formula, or in the semiclassi- 
cal expression for the spectral determinant. To extract the semiclassical spectrum 
we recall that the exact spectrum, {E n }, can be obtained from the exact counting 
function: 

oo 

N(E) = J2®(E~E n ) , (114) 

n=X 

by solving the equation 

N(E n )=n-^, n = l,2,... . (115) 

In the last equation, an arbitrarily small amount of smoothing must be applied 
to the Heavyside function. In analogy, one obtains the semiclassical spectrum 
{£-}as@: 

N BC (E%)=n-±, n = l,2,... , (116) 

where N sc is a semiclassical approximation of N. Note that N sc with which we 
start is not necessarily a sharp counting function. However, once {E^} is known, 
we can "rectify" the smooth N sc into the sharp counting function N# |J: 

oo 

N#(E)=J2 e ( E - E n)- ( 117 ) 

n=l 

The most obvious choice for N sc is the Gutzwiller trace formula truncated at 
the Heisenberg time, which is what we shall use. Alternatively, one can start 
from the regularized Berry-Keating Zeta function ( SC (E) HSfiH , and define N sc = 
N - (l/vr) Im log ( SC (E + «0), which yields N sc = N*. 

Next, in order to define a quantitative measure of the semiclassical error, 
one should establish a one-to-one correspondence between the quantal and the 
semiclassical levels, namely, one should identify the semiclassical counterparts of 
the exact quantum levels. In classically chaotic systems, for which the Gutzwiller 
trace formula is applicable, the only constant of the motion is the energy. This is 
translated into a single "good" quantum number in the quantum spectrum, which 
is the ordinal number of the levels when ordered by their magnitude. Thus, the 
only correspondence which can be established between the exact spectrum {E n } 
and its semiclassical approximation, {E^}, is 

E n <— > E s n c . (118) 
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This is to be contrasted with integrable systems, where it is appropriate to com- 
pare the exact and approximate levels which have the same quantum numbers. 

The natural scale on which the accuracy of semiclassical energy levels should 
be measured is the mean level spacing (^(i?))" 1 . We shall be interested here in the 
mean semiclassical error, and proper measures are the mean absolute difference: 



e^{E) = {d{E n )\E n -E^\ ) E (119) 



or the variance: 



eV\E) = ( (d(E n )(E n -E s n c )Y) E , (120) 

where (•) denotes averaging over a spectral interval AE centered at E. The 
interval AE is large enough so that the mean number of levels AE ■ d(E) 3> 1. 
Yet, AE is small enough on the classical scale, such that d(E) « constant over 
the interval considered. 

We shall now compare two different estimates for the semiclassical error. The 
first one is the traditional estimate: 

traditional = Q^-d^ _^ f ^OUSt , d = 2 as ft (121) 

I oo ; d ^_ o 

(c.f. section [[]). It claims that the semiclassical approximation is (marginally) 
accurate in two dimensions, but it fails to predict accurate energy levels for three 
dimensions or more. We emphasize that the traditional estimate is a qualita- 
tive error measure, emerging from global error estimate of the time propagator. 
Hence, it cannot be directly connected to either or e^ 2 -*. We mention it here 
since it is the one usually quoted in the literature. 

One may get a different estimate of the semiclassical error, if the Gutzwiller 
Trace Formula (GTF) is used as a starting point. Suppose that we have calculated 
N sc to a certain degree of precision, and we compute from it the semiclassical 



energies Eff using (|116| ). Denote by AN SC the higher order terms which were 
neglected in the calculation of N sc . The expected error in E^ can be estimated 
by including AiV sc and calculating the energy differences 5 n . That is, we consider: 

N SC (E™ + S n ) + AN sc (E s n c + 5 n )=n- 1 -. (122) 



Combining (|116|) and ( |122j ) we get (to first order in S n ): 



n ~ dN K (E«)/dE ~ d{E%) ' 1 1 

In the above we assumed that the fluctuations of N sc around its average are not 
very large. Thus, 

e«' GTF « d(E s n c )\6 n \ « AN sc (E s n c ) . (124) 
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Let us apply the above formula and consider the case in which we take for N sc 
its mean part N, and that we include in N terms of order up to (and including) 
h~ m , m < d. For AiV sc we use both the leading correction to N and the leading 
order periodic orbit sum which is formally (termwise) of order h°. Hence, 

e (l),GTF = o^-m+l) + = q ^ m in(-m+l,0)J _ ( 125 ) 

We conclude, that approximating the energies only by the mean counting function 
N up to (and not including) the constant term, is already sufficient to obtain 
semiclassical energies which are accurate to O(H ) = 0(1) with respect to the 
mean density of states. Note again, that no periodic orbit contributions were 
included in N sc . Including less terms in N will lead to a diverging semiclassical 
error, while more terms will be masked by the periodic orbit (oscillatory) term. 
One can do even better if one includes in N sc the smooth terms up to and including 
the constant term (O(h )) together with the leading order periodic orbit sum 
which is formally also O(h ). The semiclassical error is then: 



€ 



(1),GTF 



po 



Oih 1 ) . (126) 



That is, the semiclassical energies measured in units of the mean level spacing are 
asymptotically accurate independently of the dimension! This estimate grossly 



contradicts the traditional estimate ( 121|) and calls for an explanation. 



The first point that should be noted is that the order of magnitude (power of 
h) of the periodic orbit sum, which we considered above to be O(h ), is only a 
formal one. Indeed, each term which is due to a single periodic orbit is of order 
O(h ). However the periodic orbit sum absolutely diverges, and at best it is only 
conditionally convergent. To give it a numerical meaning, the periodic orbit sum 
must therefore be regularized. This is effectively achieved by truncating the trace 
formula or the corresponding spectral £ function |^6|, |84], ^9], |9(|. However, the 



cutoff itself depends on h. One can conclude, that the simple-minded estimate 
( |126|) given above is at best a lower bound, and the error introduced by the 



periodic orbit sum must be re-evaluated with more care. This point will be dealt 
with in great detail in the sequel, and we shall eventually develop a meaningful 
framework for evaluating the magnitude of the periodic orbit sum. 

The disparity between the traditional estimate of the semiclassical error and 
the one based on the trace formula can be further illustrated by the following 
argument. The periodic orbit formula is derived from the semiclassical propagator 
K sc using further approximations pf. Therefore one wonders, how can it be that 
further approximations of K sc actually reduce the semiclassical error from (VH.) 



to ( |126| )? The puzzle is resolved if we recall, that in order to obtain e^o^ 
above we separated the density of states into a smooth part and an oscillating 
part, and we required that the smooth part is accurate enough. To achieve this, 
we have to go beyond the leading Weyl's term and to use specialized methods to 
calculate the smooth density of states beyond the leading order. These methods 
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are mostly developed for billiards |82| , |91| , |92f . In any case, to obtain ep we 
have added additional information which goes beyond the leading semiclassical 
approximation. 

A direct check of the accuracy of the semiclassical spectrum using the error 
measures e- l \ is exceedingly difficult due to the exponentially large number 
of periodic orbits needed. The few cases where such tests were carried out involve 
2D systems and it was possible to check only the lowest (less than a hundred) 



levels (e.g. [60, |73|). The good agreement between the exact and the semiclassical 



values confirmed the expectation that in 2D the semiclassical error is small. In 



3D, the topological entropy is typically much larger and the direct test 

of the semiclassical spectrum becomes prohibitive. 

Facing with this grim reality, we have to introduce alternative error mea- 
sures which yield the desired information, but which are more appropriate for a 
practical calculation. We construct the measure: 

6^(E) = (\N(E)-N*(E)\ 2 ) e . (127) 

As before, the triangular brackets indicate averaging over an energy interval AE 
about E. We shall now show that 5^ faithfully reflects the deviations between the 
spectra, and is closely related to and e^ 2 \ Note, that the following arguments 
are purely statistical and apply to every pair of staircase functions. 

Suppose first, that all the differences E^ — E n are smaller than the mean 
spacing. Then, \N — N#\ is either or 1 in most of the cases (see figure |32| ). 
Hence, \N — N#\ = \N — N#\ 2 along most of the E axis. Consequently, 

5 {2) (E) w (\N(E) - N*(E)\) E , for small deviations . (128) 

The right hand side of the above equation (the fraction of non-zero contributions) 
equals e* 1 ). Thus, 

ps , for small deviations . (129) 

If, on the other hand, deviations are much larger than one mean spacing, the 
typical horizontal distance d\E— E n \ should be comparable to the vertical distance 
\N — iV#|, and hence, in this limit 

5^ ps e*- 2 ** , for large deviations . (130) 

Therefore, we expect 5^ to interpolate between e^ 1 ) and throughout the entire 
range of deviations. This behavior was indeed observed in a numerical tests which 



were performed to check the above expectations Kl. Moreover, it was shown in 



53|] that 5^ is completely equivalent to when the spectral counting functions 
are replaced with their smooth counterparts, provided that the smoothing width 
is of the order of 1 mean level spacing and the same smoothing is applied to both 
counting functions. That is, 

Smooth ~ e(2) > for a11 deviations . (131) 
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Figure 32: Illustration of \N(E) — N*(E) \ for small deviations between quantum 
and semiclassical energies: e*- 1 -* <C d^ 1 = A. The quantum staircase N(E) is 
denoted by the full line and the semiclassical staircase N#(E) is denoted by the 
dashed line. The difference is shaded. 
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In testing the semiclassical accuracy, this kind of smoothing is essential and will be 
introduced by truncating the trace formula at the Heisenberg time tn = hd. These 
properties of the measure 5^ 2 \ and its complete equivalence to for smooth 
counting functions, renders it a most appropriate measure of the semiclassical 
error. 

We now turn to the practical evaluation of To perform the energy averag- 
ing, we choose a positive window function w(E' — E) which has a width AE near 
E and is normalized by dE'w 2 (E') = 1. It falls off sufficiently rapidly so 
that all the expressions which follow are well behaved. We consider the following 
counting functions that have an effective support on an interval of size AE about 
E: 

N(E';E) = w(E'-E)N(E') (132) 
N*(E';E) ee w(E'-E)N*(E>). (133) 

The functions N and iV* are sharp staircases, since the multiplication with w 
preserves the sharpness of the stairs (it is not a convolution!). We now explicitly 
construct S^i^E) as: 

/+oo 2 
dE' N(E';E)-N*(E';E) 
■oo 

= / dE'\N(E')-N*{E')\w 2 {E'-E). (134) 

J — oo 

To obtain 6^ ooth we need to smooth N and N# over a scale of order of one mean 
spacing. One can, e.g., replace the sharp stairs by error functions. As for N#, we 
prefer to simply replace it with the original N sc , which we assume to be smooth 
over one mean spacing. That is, we suppose that N sc contains periodic orbits up 
to Heisenberg time. Hence, 

/+oo 
dE' \N smooth (E') - N sc (E')\ 2 w 2 (E' - E). (135) 
-oo 



A comment is in order here. Strictly speaking, to satisfy ( |131| ) we need to apply 
the same smoothing to A" and to AT*, and in general j\f#' smooth ^ ]V SC , but there 
are differences of order 1 between the two functions. However, since our goal 
is to determine whether the semiclassical error remains finite or diverges in the 
semiclassical limit H — > 0, we disregard such inaccuracies of order 1. If a more 
accurate error measure is needed, then more care should be practised in this and 
in the following steps. 



Applying Parseval's theorem to (|135|) we get: 

<fiooth(£) = I / df D(t; E) - D sc (t; E) (136) 
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where 

i r+oo 

D(t; E) = —= \ dE' N smooth {E';E)exp{iE't/h) (137) 

V27T J^oo 
i r+oo 

D sc (t;E) = —= / dE' N BC {E';E)exp(iE't/H). (138) 

We shall refer to D, D sc as the (regularized) quantal and semiclassical time spec- 
tra, respectively. These functions are the analogs of the length spectra D(l;k) 
used in section [5] for the billiard problem. The analogue becomes clear by in- 
voking the Gutzwiller trace formula and expressing the semiclassical counting 
function as a mean part plus a sum over periodic orbits. We have: 

N SC (E) = N(E) + ^T§- MSjW/h - ¥ /2] , (139) 

po 3\ J 



where Aj = Tj/(-nhrj^\ det(7 — Mj)\) is the semiclassical amplitude of the j'th 
periodic orbit, and Tj, Sj,Uj, Mj,rj are its period, action, Maslov index, mon- 
odromy matrix and repetition index, respectively. Then, the corresponding time 
spectrum reads: 

D sc (t;E) « D(t; E) (140) 
+ i E h 4j§l {e m)[Et+S * m ^lt + TAE)]/h)- 



TAE) 

po J v 



w 



{[t-TjWyn)}. 



In the above, the Fourier transform of w is denoted by w. It is a localized 
function of t whose width is At w h/AE. The sum over the periodic orbits 
in D sc therefore produces sharp peaks centered at times that correspond to the 
periods Tj. The term D corresponds to the smooth part and is sharply peaked 



near t = 0. To obtain ( |14Q|) we expanded the actions near E to first order: 
Sj(E') ~ Sj(E) + (E' — E)Tj(E). We note in passing, that this approximate 
expansion of Sj can be avoided altogether if one performs the Fourier transform 
over h~ l rather than over the energy. This way, an action spectrum will emerge, 
but also here the action resolution will be finite, because the range of h~ l should 
be limited to the range where d(E; h) is approximately constant. It turns out 
therefore, that the two approaches are essentially equivalent, and for billiards 
they are identical. 

The manipulations done thus far were purely formal, and did not manifestly 
circumvent the difficult task of evaluating 5g mooth . However, the introduction 



smooth" 

of the time spectra and the formula (|136|) put us in a better position than the 



original expression ( 134 ). The advantages of using the time spectra in the present 



context are the following: 
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• The semiclassical time spectrum D sc (t;E) is absolutely convergent for all 
times (as long as the window function w is well behaved, e.g. it is a Gaus- 
sian). This statement is correct even if the sum ( |140| ) extends over the 
entire set of periodic orbits! This is in contrast with the trace formula ex- 
pression for N sc (and therefore N sc ) which is absolutely divergent if all of 
the periodic orbits are included. 

• Time scale separation: As we noted above, the time spectrum is peaked 
at times that correspond to periods of the classical periodic orbits. This 
allows us to distinguish between various qualitatively different types of con- 
tributions to <5 s ( 2 oth- 

We shall now pursue the separation of the time scales in detail. We first 
note, that due to N, N sc being real, there is a t <-> (— t) symmetry in (|136|), and 



therefore the time integration can be restricted to the limits to +00: 6^ 

PCX 

la 



{1/h) ■ ■ ■ = (2/h) J °° • • •. We now divide the time axis into four intervals 



< t < At: The shortest time scale in our problem is At = h/AE. The 
contributions to this time interval are due to the differences between the 
exact and the semiclassical mean densities of states. This is an important 
observation, since it allows us to distinguish between the two sources of 
semiclassical error — the error that emerges from the mean densities and the 
error that originates from the fluctuating part (periodic orbits). Since we 
are interested only in the semiclassical error that results from the fluctuating 
part of the spectral density, we shall ignore this regime in the following. 

At < t < t erg : This is the non-universal regime |65fl , in which periodic orbits 



are still sparse, and cannot be characterized statistically. The "ergodic" 
time scale t erg is purely classical and is independent of h. 

terg < t < tn: In this time regime periodic orbits are already in the universal 
regime and are dense enough to justify a statistical approach to their pro- 
liferation and stability. The upper limit of this interval is the Heisenberg 
time £h = hd(E), which is the time that is needed to resolve the quan- 
tum (discrete) nature of a wavepacket with energy concentrated near E. 
The Heisenberg time is "quantal" in the sense that it is dependent of h: 

t u = 0{h l - d ). 

< t < 00: This is the regime of "long" orbits which is effectively truncated 
from the integration as a result of the introduction of a smoothing of the 
quantal and semiclassical counting functions, with a smoothing scale of the 
order of a mean level spacing. 

Dividing the integral ( |136| ) according to the above time intervals, we can rewrite 
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*(2) . 
smooth - 



c(2) 
smooth 



(E) 



= 5. 




D{t;E)-D sc {t; E) 



short 



long 



(141) 



(2) 

As explained above, can be ignored due to smoothing on the scale of a mean 

level spacing. The integral 5 short is to be neglected for the following reason. The 
integral extends over a time interval which is finite and independent of h, and 
therefore it contains a fixed number of periodic orbits contributions. The semi- 
classical approximation provides, for each individual contribution, the leading 
order in h, and therefore |8(| we should expect: 



6. 



(2) 
short 



as fi 



0. 



(142) 



Our purpose is to check whether the semiclassical error is finite or divergent as 
h — > 0, and to study whether the rate of divergence depends on dimensionality. 
Equation ( 142|) implies that 5^ ort cannot affect 5^ in the semiclassical limit and 
we shall neglect it in the following. 
We remain with: 



<5 (2) « 

smooth 



§ (2) 
u m > 



(143) 



which will be our object of interest from now on. 

The fact that tu is extremely large on the classical scale renders the calculation 
of all the periodic orbits with periods less than an impossible task. However, 
sums over periodic orbits when the period is longer than t erg tend to meaningful 

(2) 

limits, and hence, we would like to recast the expression for S m in the following 
way. Write 5m as: 



dt 



dt 



b{t) - D sc (t) 



bit) 



D(t) - b sc {t) 



b(t) 



2 f til 



h 



dt 



'ere 



D(t) ) x C(t) 



(144) 



(145) 



(146) 



/•tn 

/ envelope x correlation 

J tore: 



where the parametric dependence on E was omitted for brevity. The smoothing 
over t is explicitly indicated to emphasize that one may use a statistical inter- 
pretation for the terms of the integrand. This is so because in this domain, the 
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density of periodic orbits is so large, that within a time interval of width h/AE 
there are exponentially many orbits whose contributions are averaged due to the 
finite resolution. 

We note now that we can use the following relation between the time spectrum 
and the spectral form factor K(t): 

Kir) 

/* d t = ffUr (147) 

where r = t/tn is the scaled time. The above form factor is smoothed according 
to the window function w. Hence: 

e^JL/v™. (148 ) 

** 7"erg 

For generic chaotic systems we expect that K(t) agrees with the results of RMT 
in the universal regime r > r erg |24|, H |65|1 . Therefore 

t < K{t) < gr for r crg < r < 1 , (149) 

where = 1 for systems which violate time reversal symmetry, and g = 2 if 

c(2) 

time reversal symmetry is respected. This implies that the evaluation of <5 smooth 
reduces to 

tfO „ 9 f 1 ,_C(r) 



Coot, -2^/ dr^, (150) 

where we took the upper bound gr for K{t). The dependence on h in this 
expression comes from the lower integration limit which is proportional to Ti d ~ l 
as well as from the implicit dependence of the function C on h. 

Formula ( |150j ) is our main theoretical result. However, we do not know how 
to evaluate the correlation function C (r) from first principles. The knowledge of 
the h corrections to each of the terms in the semiclassical time spectrum is not 
sufficient since the resulting series which ought to be summed is not absolutely 
convergent. Therefore we have to recourse to a numerical analysis, which will 
be described in the next section. The numerical approach requires one further 
approximation, which is imposed by the fact that the number of periodic orbits 
with t < tft is prohibitively large. We had to limit the data base of periodic orbits 
to the domain t < t cpn with t erg £ cpu tn- The time t cpu has no physical origin, 
it represents only the limits of our computational resources. Using the available 
numerical data we were able to compute C(t) numerically for all t eTg < t < t cpu 
and we then extrapolated it to the entire domain of interest. We consider this 
extrapolation procedure to be the main source of uncertainty. However, since the 
extrapolation is carried out in the universal regime, it should be valid if there are 
no other time scales between t eTg and tn. 



99 



6.2 Numerical results 



We used the formalism and definitions presented above to check the accuracy of 
the semiclassical spectra of the 2D and 3D Sinai billiards. The most important 
ingredient in this numerical study is that we could apply the same analysis to 
the two systems, and by comparing them to give a reliable answer to the main 
question posed in this work, namely, how does the semiclassical accuracy depend 
on dimensionality. 

The classical dynamics in billiards depends on the energy (velocity) trivially, 
and therefore the relevant parameter is the length rather than the period of the 
periodic orbits. Likewise, the quantum wavenumbers k n are the relevant variables 
in the quantum description. From now on we shall use the variables (/, k) instead 
of {t,E), and use "length spectra" rather than "time spectra". The semiclassical 
limit is obtained for k — ► oo and O(h) is equivalent to 0{k~ 1 ). Note also that 
for a billiard N(k) ~ Ak d where A is proportional to the billiard's volume. 

We start with the 2D Sinai billiard, which is the free space between a square 
of edge S and an inscribed disc of radius R, with 2R < S. Specifically, we use 
S = 1 and R = 0.25 and consider the quarter desymmetrized billiard with Dirich- 
let boundary conditions for the quantum calculations. The quantal database con- 
sists of the lowest 27645 eigenvalues in the range < k < 1320, with eigenstates 
which are either symmetric or antisymmetric with respect to reflection on the 
main diagonal. The classical database consists of the shortest 20273 periodic or- 
bits (including time reversal, reflection symmetries and repetitions) in the length 
range < / < 5. For each orbit, the length, the stability determinant and the 
reflection phase were recorded. The numerical work is based on the quantum 
spectra and on the classical periodic orbits which were computed by Schanz and 
Smilansky [0, g for the 2D billiard. 

We begin the numerical analysis by demonstrating numerically the correctness 
of equation ( |142[ ). That is, that for each individual contribution of a periodic 
orbit, the semiclassical error indeed vanishes in the semiclassical limit. In figure 



53] we plot \D — D sc \ for I = 0.5 as a function of k. This length corresponds to the 
shortest periodic orbit, that is, the one that runs along the edge that connects the 
circle with the outer square. For D sc we used the Gutzwiller trace formula. As is 
clearly seen from the figure, the quantal-semiclassical difference indeed vanishes 
(approximately as A; -1 ), in accordance with (|142|) . We emphasize again, that 
this behavior does not imply that 8^ 2 > vanishes in the semiclassical limit, since 
the number of periodic orbits included depends on k. It implies only that 5^ 0Tt 
vanishes in the limit, since it consists of a fixed and finite number of periodic orbit 
contributions. We should also comment that penumbra corrections to individual 
grazing orbits introduce errors which are of order k' 1 with < 7 < 1 
However, since the definition of "grazing" is in itself k dependent, one can safely 
neglect penumbra corrections in estimating the large k behavior of S^ OTt . 

We now turn to the main body of the analysis, which is the evaluation of 
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k 

Figure 33: The absolute difference between the quantal and the semiclassical 
(Gutzwiller) length spectra for the 2D Sinai billiard at I = 0.5. This length 
corresponds to the shortest unstable periodic orbit. The average log-log slope is 
about —1.1, indicating approximately k~ x decay. The data were averaged with a 
Gaussian window. 
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(2) 

dm for the 2D Sinai billiard. Based on the available data sets, we plot in figure 



34] the function C(l; k) in the interval 2.5 < I < 5 for various values of k. One 
can observe, that as a function of I the functions C(l; k) fluctuate in the interval 
for which numerical data were available, without exhibiting any systematic mean 
trend to increase or to decrease. We therefore approximate C(l; k) by 

C(l;k) ss const ■ f(k) = C &vg {k). (151) 

As mentioned above, we extrapolate this formula in I up to the Heisenberg length 
L H = 27rJ(/c) and using (|150|) we obtain: 



*££l = %^ M WL erg ) = C avg (k) 0(ln k). (152) 

The last equality is due to L H = Oik^ 1 ). To evaluate C avg (/c) we averaged 
C(l; k) over the interval L eTg = 3.5 < / < 5 = L cpu and the results are shown in 
figure [35| We choose L erg = 3.5 because the density of periodic orbits is large 
enough for this length (see figure |4j) to expect universal behavior of the periodic 
orbits. (For the Sinai billiard described by flow the approach to the invariant 
measure is algebraic rather than exponential f|0, |39|, and thus one cannot have 



a well-defined L evg . An any rate, the specific choice of L erg did not affect the 
results in any appreciable way.) Inspecting C avg (/c), it is difficult to arrive at firm 
conclusions, since it seems to fluctuate around a constant value up to k ss 900 and 
then to decline. If we approximate C avg (k) by a constant, we get a "pessimistic" 
value of 5^: 

8^'™ h (k) = 0(\nk) = 0(\nh) "pessimistic" (153) 

while if we assume that C avg (A;) decays as a power-law, C avg (k) = > 0, 

then 

€Sh(*0 = 0{k-P\nk) — "optimistic" . (154) 
Collecting the two bounds we get: 

O(k^lnk) < 5^Z(k) < OQnk) . (155) 

Our estimates for the 2D Sinai billiard can be summarized by stating that the 
semiclassical error diverges no worse than logarithmically (meaning, very mildly) . 
It may well be true that the semiclassical error is constant or even vanishes in the 
semiclassical limit. To reach a conclusive answer one should invest exponentially 
larger amount of numerical work. 

There are a few comments in order here. Firstly, the quarter desymmetrization 
of the 2D Sinai billiard does not exhaust its symmetry group, and in fact, a 
reflection symmetry around the diagonal of the square remains. This means, 
that the spectrum of the quarter 2D Sinai billiard is composed of two independent 
spectra, which differ by their parity with respect to the diagonal. If we assume 
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Figure 34: The functions C(l; k) for quarter 2D Sinai billiard S = 1,R = 0.25 
with Dirichlet boundary conditions. The window w(k' — k) was taken to be a 
Gaussian with standard deviation a = 60. We averaged C(l; k) over /-intervals of 
« 0.2 in accordance with ( |145| ) to avoid sharp peaks due to small denominators. 
The averaging, however, is fine enough not to wash out all of the features of 
C(l;k). The vertical bars indicate the locations of primitive periodic orbits, and 
the daggers indicate the locations of the bouncing-ball families. 
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that the semiclassical deviations of the two spectra are not correlated, the above 
measure is the sum of the two independent measures. It is plausible to assume also 
that both spectra have roughly the same semiclassical deviation, and thus Smooth 
is twice the semiclassical deviation of each of the spectra. Secondly, we recall 
that the 2D Sinai billiard contains "bouncing-ball" families of neutrally stable 



periodic orbits |T5|, |37|, We have subtracted their leading-order contribution 
from D such that it includes (to leading order) only contributions from generic, 
isolated and unstable periodic orbits. This is done since we would like to deduce 
from the 2D Sinai billiard on the 2D generic case in which the bouncing-balls 
are not present. (In the Sinai billiard, which is concave, there are also diffraction 



effects |81j, |3§H , but we did not treat them here.) Thirdly, the analogue of (|147f) 
for billiards reads: 

D(l) 2 ) dl = d£ (156) 



4tt 2 £ 2 



when £ = l/Lu- In figure |36| we demonstrate the compliance of the form factor 
with RMT GOE using the integrated version of the above relation, and taking 
into account the presence of two independent spectra. Finally, it is interesting to 
know the actual numerical values of 5 s ^ h (k) for the k values that we considered. 



We carried out the computation, and the results are presented in figure 37. One 
observes that for the entire range we have ^sniooth(^) ~ 0-1 ^ ^ which is very 
encouraging from an "engineering" point of view. 

We now turn to the analysis of the 3D Sinai billiard. We use the longest 
quantal spectrum (R = 0.2, Dirichlet) and the classical periodic orbits with 
length < I < 5. 

To treat the 3D Sinai billiard we have to somewhat modify the formalism 
which was presented above. This is due to the fact that in the 3D case the 
contributions of the various non-generic bouncing-ball manifolds overwhelm the 
spectrum [^, and unlike the 2D case, it is difficult to explicitly eliminate 
their (leading-order) contributions (c.f. the discussion in section |5.4| ). Since our 
goal is to give an indication of the semiclassical error in generic systems, it is 
imperative to avoid this dominant and non- generic effect. 

We shall use the mixed boundary conditions, which were discussed in section 



5~q and were shown to largely filter the bouncing-ball effects. Specifically, we 
consider d (c.f. ( |113| )) for our purposes. Let us construct the weighted counting 
function: 

N(k) = / dk' d(k') = V v n Q(k - k n ) . (157) 

The function N is a staircase with stairs of variable height v n . As explained 
above, its advantage over N is that it is semiclassically free of the bouncing 



balls (to leading order) and corresponds only to the generic periodic orbits [[33 
Similarly, we construct from d sc the function N sc . Having defined N, N sc , we 
proceed in analogy to the Dirichlet case. We form from N, N sc the functions 
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Figure 36: Verification of equation (|156|) for the quarter 2D Sinai billiard. We 
plot /(£) = d£' K 2 and compare the quantum data with RMT. The 
minimal £ corresponds to L erg = 3.5. The integration is done for smoothing, 
and we fix the upper limit to avoid biases due to non-universal regime. Note the 
logarithmic scale. 
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Figure 37: The numerical values of Smooth f° r the quarter 2D Sinai billiard. We 

(2) 

included also the contribution S s ^ ort of the non-universal regime. The contribu- 

(2} (2) 

tions from the time interval t erg < t < t cpn are contained in <5m,c P u, and 5^ ext is 
the extrapolated value for t cpu < t < tn (refer to equation (|141| ) and to the end 
of subsection [67l] ) . 
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N, N sc , respectively, by multiplication with a window function w(k' — k) and 
then construct the measure 5^ as in (|134j) . The only difference is that the 
normalization of w must be modified to account for the "velocities" v n such as: 

d-\k)^v 2 n \w(k n - k)\ 2 = 1. (158) 

n 

The above considerations are meaningful provided the "velocities" v n are nar- 
rowly distributed around a well-defined mean v(k), and we consider a small 
enough fc-interval, such that v(k) does not change appreciably within this in- 
terval. Otherwise, 5^ is greatly affected by the fluctuations of v n (which is 
undesired) and the meaning of the normalization is questionable. We shall check 
this point numerically. 

To demonstrate the utility of the above construction using the mixed bound- 
ary conditions, we return to the 2D case. We set k = 1007T, and note that 
the spectrum at our disposal for the mixed case was confined to the interval 
< k < 600. First, we want to examine the width of the distribution of the 



u n 's. In figure [38| we plot the ratio of the standard deviation of v n to the mean, 
averaged over the /c-axis using a Gaussian window. We use the same window also 
in the calculations below. The observation is that the distribution of v n is mod- 
erately narrow and the width decreases algebraically as k increases. This justifies 
the use of the mixed boundary conditions as was discussed above. One also needs 
to check the validity of (|156|) , and indeed we found compliance with GOE also 
for the mixed case (results not shown). We next compare the functions C(l; k) 
for both the Dirichlet and the mixed boundary conditions. It turns out, that also 
in the mixed case the functions C(l;k) fluctuate in / with no special tendency 
(not shown). The averages C avg (A;) for the Dirichlet and mixed cases are com- 



pared in figure |39[ The values in the mixed case are systematically smaller than 



in the Dirichlet case which is explained by the efficient filtering of tangent and 



close to tangent orbits that are vulnerable to large diffraction corrections |81 
However, from k = 250 on the two graphs show the same trends, and the values 
of C avg in both cases are of the same magnitude. Thus, the qualitative behav- 
ior of dgmooth * s shown to be equivalent in the Dirichlet and mixed cases, which 

(2) 

gives us confidence in using <5 sn 4, oth together with the mixed boundary conditions 
procedure. 

We finally applied the mixed boundary conditions procedure to compute 
Smooth f° r the desymmetrized 3D Sinai with S = 1,R = 0.2 and set k = 100. 
We first verified that also in the 3D case the velocities v n have a narrow distri- 
bution — see figure Next, we examined equation (|156|) using quantal data, 



and discovered that there are deviations form GOE (figure fH]). We have yet 
no satisfactory explanation of these deviations, but we suspect that they are 
caused because the ergodic limit is not yet reached for the length regime under 
consideration due to the effects of the infinite horizon which are more acute in 
3D. Nevertheless, from observing the figure as well as suggested by semiclassical 
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Figure 38: Calculation of Q = a/ (v%) — ( v n) 2 /\(v n )\ for quarter 2D Sinai billiard 
(up) and for the desymmetrized 3D Sinai billiard (down). 



109 




k 

Figure 39: Comparison of C avg (/c) for Dirichlet and mixed boundary conditions 
for the quarter 2D Sinai billiard. We used a Gaussian window with a = 40. 
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arguments, it is plausible to assume that K(£) oc £ for small £. Hence, this devi- 
ation should not have any qualitative effect on 5^ according to (|150|) . Similarly 
to the 2D case, the behavior of the function C(l; k) is ffuctuative in I, with no 
special tendency (figure f|I]). If we average C(l;k) over the universal interval 
L erg = 2.5 < / < L cpu = 5, we obtain C avg (k) which is shown in figure ^2|. The 
averages C avg (fc) are fluctuating with a mild decrease in k, and therefore we can 
again conclude that 

0(k-?\nk)<5i±™ h <0(lnk) (159) 

where the "optimistic" measure (leftmost term) corresponds to C avg (fc) = Oik^ 13 ), (3 > 
0, and the "pessimistic" one (rightmost term) is due to C avg (k) = const. In other 
words, the error estimates ( |155| , |159| ) for the 2D and the 3D cases, respectively, 
are the same, and in sharp contrast to the traditional error estimate which pre- 
dicts that the errors should be different by a factor 0(h~ l ). On the basis of our 
numerical data, and in spite of the uncertainties which were clearly delineated, 
we can safely rule out the traditional error estimate. 

Our main finding is that the upper bound on the semiclassical error is a 
logarithmic divergence, both for a generic 2D and 3D systems (equations ( |155|) , 
( |159|) ). In this respect, there are a few points which deserve discussion. 

c(2) 

To begin, we shall try to evaluate 5 smoo th usm g the explicit expressions for the 
leading corrections to the semiclassical counting function of a 2D generic billiard 
system, as derived by Alonso and Gaspard [0: 



A; 



N(k) = N(k) + — sin 



kL 3 + QjL + 0{l/k 2 ) 



(160) 



where Aj are the standard semiclassical amplitudes, Lj are the lengths of periodic 
orbits and Qj are the fc-independent amplitudes of the 1/k corrections. The Q/s 
are explicitly given in 0. We ignored in the above equation the case of odd 
Maslov indices. If we calculate from N(k) the corresponding length spectrum 
D(l; k) using a (normalized) Gaussian window w(k' — k) — (1/ ^/na 2 ) exp[— {k' — 
k) 2 /(2a 2 )], we obtain: 



D(l;k) 



2^F 



E 



L 



{l-Lj) 2 ^ _ ik{l+Lj)+i- 



k e 



j u 



(161) 

In the above we regarded the phase e l< ^ i//fc as slowly varying. The results of Alonso 
and Gaspard [0 suggest that the Qj are approximately proportional to the length 
of the corresponding periodic orbits: 



Qj ~ QLj. 

We can therefore well-approximate D as: 

l _V°_ „-iQl/k A o r 
2^F 6 ^ Lj [ - " 



D(l;k) 



iQl/k -q 



sc-GTF 



(162) 



(163) 
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Figure 40: Check of equation (|156|) for the desymmetrized 3D Sinai billiard. The 



minimal £ corresponds to L erg = 2.5. The function /(£) is defined as in figure |36 
Note the logarithmic scale. 
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Figure 41: The functions C(l] k) for desymmetrized 3D Sinai billiard S — 1,R — 
0.2 with mixed boundary conditions. We took a Gaussian window with a = 20, 
and smoothed over /-intervals of w 0.3. The upper vertical bars indicate the 
locations of primitive periodic orbits. 
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Figure 42: Averaging in / of C(l; k) for 3D Sinai billiard as a function of k. The 
averaging was performed in the interval L erg = 2.5 < / < 5 = L cpn . 
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where Z) sc _gtf is the length spectrum which corresponds to the semiclassical 
Gutzwiller trace formula for the counting function (without 1/k corrections). We 
are now in a position to evaluate the semiclassical error, indeed: 

*£Lh(*) = 2 dl D(l;k)-D sc ^ GTF (l-kf 



/ §)| J& ( i ;*)f- ( 164 ) 



dl sin 2 



If we use equation (|156|) and K(l) ~ gl/Ln (which is valid for I < Lh for chaotic 
systems), we get: 

c(2) ^ „ 2 9 [ Ln dl ^ 2 (Ql\ _ 2g ^) ^sin 2 (t) 



^k)»%r***^. (167) 



€U W « J / y sin 2 =r = J / dt— ^ . (165) 

^ A min f tH J QLm .J(2k) t 

For fc^oowe have that 

/ dt^IZ^ / dt ■ t = 0(l/k 2 ) (166) 

Jo * Jo 

which is negligible, hence we can replace the lower limit in ( |165| ) with 0: 

vr 2 Jo t 

(2) 

This is the desired expression. The dimensionality enter in <5g mooth (fc) only through 
the power of k in Lh. 

Let us apply equation (|167|) to the 2D and the 3D cases. For 2D we have to 
leading order that Lh = Ak, where A is the billiard's area, thus, 

<££-(*) - | r\t S ^ = const = 0(f) (168) 

which means, that the semiclassical error in 2D billiards is of the order of the mean 
spacing, and therefore the semiclassical trace formula is (marginally) accurate and 
meaningful. This is compatible with our numerical findings. 

For 3D, the coefficients Qj were not obtained explicitly, but we shall assume 
that they are still proportional to Lj (equation ( |162|) ) and therefore that ( |167| ) 
holds. For 3D billiards Lh = (V/7r)k 2 to leading order, where V is the billiard's 
volume. Thus the upper limit in ( |167| ) is QVk/(27r) which is large in the semi- 
classical limit. In this case, we can replace sin 2 (t) with its mean value 1/2 and 
the integrand becomes essentially 1/t which results in: 

*2$c-(*)=O(lnA0. (169) 
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That is, in contrast to the 2D case, the semiclassical error diverges logarithmically 
and the semiclassical trace formula becomes meaningless as far as the prediction 
of individual levels is concerned. This statement is compatible with our numer- 
ical results within the numerical dispersion. However, it relies heavily on the 
assumption that Qj ~ QLj, for which we can offer no justification. We note in 
passing, that the logarithmic divergence persists also for d > 3. 

Another interesting point relates to integrable systems. It can happen, that for 
an integrable system it is either difficult or impossible to express the Hamiltonian 
as an explicit function of the action variables. In that case, we cannot assign to the 
levels other quantum numbers than their ordinal number, and the semiclassical 
error can be estimated using S^. However, since for integrable systems K(r) = 1, 
we get that: 

,(2),int dr £fcl ( 17 Q) 

°smooth ~ 2vr 2 aT T 2 ■ ^ llU ) 

J 7"erg 

Therefore, for deviations which are comparable to the chaotic cases, Cir) = 0(1), 
we get (^smooth = 0{h l ~ d ) which is much larger than for the chaotic case and 
diverges for d > 1. 

The formula ( |l50p for the semiclassical error contains semiclassical informa- 
tion in two respects. Obviously, C(r), which describes the difference between 
the quantal and the semiclassical length spectra, contains semiclassical informa- 
tion. But also the fact that the lower limit of the integral in (|150|) is finite is 
a consequence of semiclassical analysis. If this lower limit is replaced by 0, the 
integral diverges for finite values of h. Therefore, the fact that the integral has a 
lower cutoff, or rather, that D is exactly below the shortest period, is a crucial 
semiclassical ingredient in our analysis. 

Finally, we consider the case in which the semiclassical error is estimated with 
no periodic orbits taken into account. That is, we want to calculate (\N(E) — 
N(E)\ 2 ) e which is the number variance E 2 (x) for the large argument x = AE d(E) ^> 



erg; 



1. This implies C{r) = 1, and using ([1501) we get that <^ ooth = £/(2tt 2 ) ln(i H /£ 
which in the semiclassical limit becomes g/(27r 2 ) ln(tn) = 0(lnh). This result is 
fully consistent and compatible with previous results for the asymptotic (satura- 
tion) value of the number variance £ 2 (see for instance |95], |96|). It implies 



also that the pessimistic error bound (|153| ) is of the same magnitude as if peri- 
odic orbits were not taken into account at all. (Periodic orbits improve, however, 
quantitatively, since in all cases we obtained C avg < 1.) Thus, if we assume that 
periodic orbit contributions do not make N sc worse than N, then the pessimistic 
error bound 0(lnh) is the maximal one in any dimension d. This excludes, in 
particular, algebraic semiclassical errors, and thus refutes the traditional estimate 

o(h 2 - d ). 
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7 Semiclassical theory of spectral statistics 



In section |] we studied several quantal spectral statistics of the Sinai billiard 
and have shown that they can be reproduced to a rather high accuracy by the 
predictions of Random Matrix Theory (RMT). In the present section we would 
like to study the spectral two-point correlation function in the semiclassical ap- 
proximation, and to show how the classical sum rules and correlations of periodic 
orbits, which were defined in section £|, can be used to reconstruct, within the 
semiclassical approximation, the predictions of RMT. 

The starting point of the present discussion is the observation that the semi- 



classical spectrum can be derived from a secular equation of the form [E6], |83 



Z sc {k) = det(J - S(k)) = , (171) 

where S(k) is a (semiclassically) unitary matrix which depends parametrically 
on the wavenumber k. In the semiclassical approximation, the unitary operator 
S(k) can be considered as the quantum analogue of a classical Poincare mapping, 
which for billiard systems in d dimensions, is the classical billiard bounce map. 
The dimension N(k) of the Hilbert space on which S(k) acts, can be expressed 
within the semiclassical approximation, in terms of the phase-space volume of 
the Poincare section M. as follows: 

N(k) = [Af(k)} , J\f(k) = , (172) 

where [•] stands for the integer value. For a billiard in two dimensions N(k) = 
Ck/n, where C is the circumference of the billiard. In the case of the fully 
desymmetrized 3D Sinai billiard, for which we consider the sphere return map, 
M{k) = k 2 R 2 /48. The reason why we defined the smooth function M{k) will 
become clear in the sequel. 

The eigenvalues of S(k) are on the unit circle: {exp(i8i(k))}^}^ . If for a 
certain k, one of the eigenphases is an integer multiple of 2ir, then equation ( |171| ) 
is satisfied, and this value of k belongs to the spectrum. Because of this connection 
between the billiard spectrum on the k axis and the eigenphase spectrum on the 
unit circle, the statistics of fc-intervals can be read off the corresponding statistics 
of the eigenphase intervals averaged over an appropriate ^-interval where N(k) is 
constant [|6|, |84| . For this reason, it is enough to study the eigenphase statistics, 
and if they can be reproduced by the predictions of RMT for the relevant circular 
ensemble, the wavenumber spectral statistics will conform with the prediction of 
RMT for the corresponding Gaussian ensemble. 

The spectral density of the matrix S(k) can be written as: 

N(k) , oo 

d qm (£; k) = J2 6 ( 6 ~ = ~2 + 27 S ( e ~ mV ^ + e^tr^D • (173) 

1=1 71 71 n=l 
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The corresponding two-point correlation function is derived by computing: 

CM = % /jT ^ d qm (9 + | k) d qm (0 - |; k)\ , (174) 

where (•) denotes an average over a fc-interval where N(k) takes the constant 
value N. The two-point spectral form factor is defined as the Fourier coefficients 
of C 2 (t]), and by substituting (|173|) in (|174|) , one finds that they are equal to 



^ (\trS n (k)\ 2 ) . RMT provides an explicit expression: 

l^ n WI 2 ) R MT = ^(^y) > (175) 



N 

1 / ,„, , ■-», ... / // 

N 



where (3 is the standard ensemble label |p9| . The most important fact to be 
noticed is that n, the "topological time", is scaled by N, which plays here the 
role of the Heisenberg time. For a Poisson ensemble: 

i<|tr^(A;)| 2 > Poisson = l. (176) 

From now on we shall be concerned with the Circular Orthogonal Ensemble 
(COE: (3 = 1). We define requivn/N. The function Kcoe{t) is a monotonically 
increasing function which starts as 2r near the origin, and bends towards its 
asymptotic value 1 in the vicinity of r = 1. For an explicit expression consult, 
e.g. |7^| . Our aim is to show that the semiclassical expression for (\trS n (k)\ 2 ) 



reproduces this behaviour when the correlations of periodic orbits are properly 
taken into account. 

Recalling that the unitary matrix S(k) is the quantum analogue of the Poincare 
map, one can express trS n (k) in terms of the n-periodic orbits of the mapping. 
If the semiclassical mapping is hyperbolic, and the billiard bounce map is con- 



sidered, one gets (72 1 



trS n (k) « V ^ _ e ^(_i)^ . (177) 

~ n |det(I-M0|i 

Here V n is the set of all n-periodic orbits of the bounce map, n p j is the period of 
the primitive orbit of which the n-periodic orbit is a multiple. The monodromy 
matrix is denoted Mj, Lj is the length, and bj is the number of bounces from the 
boundaries (for a Dirichlet boundary condition). Note that when the Poincare 
section consists of a part of the boundary (as is the case for the sphere return 
map in the 3D Sinai billiard), bj can be different from n. Recalling the definition 
of the classical density d c \{l; n) ( f72|) in subsection [44], and realizing that the pre- 
exponential factors are just the Aj coefficients (0), we deduce that within the 
semiclassical approximation, 

/*27T /»00 

/ e ine d qm (6; k) d6 = trS n (k) w / e ikl d cl (l; n) dl . (178) 
Jo Jo 
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This equation is of fundamental importance, because it expresses the duality 
between the quantum mechanical spectral density and the classical length density 



via their Fourier transforms [11|. Hence, the spectral form factors of the classical 



and the quantum spectral distributions are also related by: 

2 \ 

(179) 



l<|tr^(A;)| 2 > = l(K cl (fc;n)) = l 



We have shown already in section |4.4j that the length spectrum as defined by 
the classical density (|72D contains non-trivial correlations. They appear on a 
scale A(n; R) which is inversely proportional to the value of k where the classical 
correlation function approaches it asymptotic value gn. What remains to be 
seen now is the extent by which the semiclassical expression ( |179j ) reproduces the 
expected universal scaling and the detailed functional dependence on the scaled 
topological time r = n/N as predicted by RMT. 

The large k limit of K c i(k;n) was written explicitly in ([F^) and verified nu- 
merically: 

K d (k — > oo; n) « (n p g p ) ■ U{n) « 2n . (180) 
This limit corresponds to the limit — * so that 

i<|trSW>-|, (181) 



which is identical to the behavior of Kqoe(t) in the small r limit [59]. Therefore, 
the classical uniform coverage of phase space guarantees the adherence to RMT 
in the limit r — > 0. This result was derived originally by Berry in his seminal 
paper 0. It is the "diagonal approximation" which can be used as long as the 
range of k values is larger than A(n; R)~ x . In other words, this approximation is 
valid on the scale on which the classical length spectrum looks uncorrelated. This 
observation shows that the domain of validity of the diagonal approximation has 
nothing to do with the "Ehrenfest time" , sometimes also called the "log h time" . 
Rather, it depends on the correlation length in the classical spectrum A(n; R), as 
displayed by the classical form factor. 

Given the classical correlation function, K c \(k;n), it cannot be meaningfully 
compared to the COE result at all values of the parameters. This is because 
once N(k) = 1, one cannot talk about quantum two-point correlations, since 
the spectrum consists of a single point on the unit circle. In other words, this 
is the extreme quantum limit, where the Hilbert space consists of a single state. 
Therefore, the fc-values to be used must exceed in the case of the 3D Sinai billiard 
^min = V9Q/R, which corresponds to Af(k) = 2. Hence, the values of r = n/N 
which are accessible are restricted to the range < r < n/2. 



In figure 43 we summarize our numerical results by comparing the the form 



factor obtained from periodic orbit theory K c \ with the theoretical RMT predic- 
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tion -Kcoe- What we actually show is the running average, 



C(t) = - [ T dr'-K cl (r'), (182) 
t Jo n 

where K c \{t) = K c \(k(r;n);n). The corresponding COE curve (c.f. equations 
( |175|) and ( |179|) ) is given by: 

C(r) = - [ T dr' K C oe(t') . (183) 
T Jo 

The "diagonal approximation" curve is obtained by replacing K coe (t) by 2t, 
namely, classical correlations are ignored. The data sets which were chosen are 
those for which sufficiently many periodic orbits were computed so that the sum 
rule U (n; I) ~ 1 was satisfied. We did not include the n — 1 data because they 
are non-generic. As clearly seen from the figure, the data are consistent with the 
RMT expression and they deviate appreciably from the diagonal approximation. 
This is entirely due to the presence of classical correlations, and it shows that the 
classical correlations are indeed responsible for the quantitative agreement. Note 
also that the data represents four different combinations of n and R, which shows 
that the classical scaling is indeed consistent with the universal scaling implied 
by RMT. In figure [P] we present essentially the same data, but integrated and 
plotted using the variable k, similarly to section |j. The integration started at 
k min for a meaningful comparison with RMT. Again, we observe the quantitative 
agreement, which is especially good for the higher n values (n = 3, 4). 

In section [O] we showed that the classical correlations originate to a large 
extent from the f2(W) families of periodic orbits. Moreover, the form factor which 
was calculated by neglecting cross-family contributions was much smoother than 
the original one. It is therefore appealing to take advantage of this smoothness 
and compare the numerical and theoretical form factors themselves instead of 
their running averages. We define: 

K cl (N; n) = (K cl (k; n)) N = k{N + x) _ k{N) J dtf K cl (k'; n) (184) 



(TV) 



we 



which is the semiclassical ensemble average of the form factor. In figure |45 
compare K c i(N; n) with N -KcoE,{n/N). The classical form factor included intra- 
family contribution only, and we multiplied it by a factor such that asymptotically 
it will match the theoretical value 2n. This factor compensates for the partial 
breaking of time-reversal symmetry and for the fact that the classical saturation 
is to values slightly below 2n for the n's under consideration. One observes that 
the agreement is quite good, and in any case the classical form factor is sharply 
different from the diagonal approximation, meaning that classical correlations 
are important. In figure [P] we present the same results with r = n/N as the 
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Figure 43: Comparing the classical form factor with the universal RMT predic- 
tions for various cases of the 3D SB. 
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variable. It again shows that the classical from factor agrees with the COE 
expression beyond the validity range of the diagonal approximation. The range 
of r where a good agreement is observed increases with n as expected, but the 
estimated domain of valid comparison r < 2n seems to be too optimistic. 

In summary we can say that the present results show that the semiclassical 
theory based on the Gutzwiller trace formula is capable to reproduce the COE 
form factor beyond the "diagonal approximation". To do this, one has to in- 
clude the classical correlations in the way which was done here, and once this 
is done, there is no need to augment the theory by uncontrolled "higher order" 
or "diffractive" corrections as was done in |85| , |86| and by others. The results 
obtained in the present section are corroborated by a recent analysis of periodic 
orbit lengths correlations in billiards constructed from octagonal modular do- 



mains in the hyperbolic plain [f74[ . The same quality of agreement was obtained 
between the classical form factor and the corresponding RMT result. These bil- 
liards are in two dimensions, and therefore the scaling laws depend differently on 
k, and the fact that the resulting scaled quantities agree with the expressions de- 
rived from RMT gives further support to the line of thinking developed here. We 
have grounds to believe that the classical correlations are universal in hyperbolic 
systems, and have to do with the self-similar organization of the set of periodic 
orbits. The previous numerical studies which were conducted also on different 



systems support this conjecture H, IT 
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8 Summary 



In the present paper we tried to provide a complete description of a paradigmatic 
three dimensional quantum system which is chaotic in the classical limit — the 
three dimensional Sinai billiard. This study is called for especially because most 
of the detailed investigations in the field were carried out for systems in two 
dimensions. 

Our main purpose in this study was to emphasize and clarify issues which 
are genuinely related to the three dimensional character of the system. The 
question which concerned us most was whether the semiclassical approximation 
- the main theoretical tool in the field — is sufficiently accurate for the spectral 
analysis of systems in three dimensions. 

We were able to obtain accurate and extensive data bases for the quantum en- 
ergy levels and for the classical periodic orbits. These allowed us to check various 
properties of the quantum spectrum, and in particular to study the applicabil- 
ity of the semiclassical approximation. The main conclusion from our work is 
that contrary to various expectations, the semiclassical accuracy, measured in 
units of the mean spacing, does not diverge as a h 2 ~ d . Our numerical tests and 
analytical arguments indicate an error margin which at worst diverges weakly 
(logarithmically) with ft. 

One of the main problems which we had to overcome was how to separate 
the generic features which are common to all chaotic systems, from the system 
specific attributes, which in the present case are the "bouncing ball" manifolds 
of periodic orbits. We should emphasize that in d dimensions the bouncing- 
ball manifolds contribute terms of order k^' 1 ^ 2 , which are much larger than the 
order 1 contributions due to generic periodic orbits. Hence it is clear that as the 
dimension increases, the extracting of generic features becomes more difficult, and 
one has also to control higher h corrections, such as, e.g., diffraction corrections 
to the bouncing ball contributions. We developed a method to circumvent some 
of these difficulties which was sufficient for the 3D Sinai billiard case, namely, we 
focused on the derivative of the spectrum with respect to the boundary condition. 
This method is a powerful means which can also be used in other instances, where 
non-generic effects should be excluded. 

One of the issues which are essential to the understanding of trace formulae 
and their application, was first mentioned by Gutzwiller in his book, under the 
title of the "third entropy" M. Gutzwiller noticed that in order that the series 
over periodic orbits can be summed up (in some sense) to a spectral density 
composed of 6 functions, the phases of the contributing terms should have very 
special relations. The more quantitative study of this problem started when 
Argaman et al. defined the concept of periodic orbit correlations. The dual 
nature of the quantum spectrum of energies and the classical spectrum of periodic 
orbit was further developed in [IT] . It follows that the universality of the quantum 
spectral fluctuations implies that the correlation length in the spectrum of the 
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classical actions depends on the dimensionality in a specific way. This was tested 
here for the first time, and the mechanism which induces classical correlations 
was discussed. 

Our work on the Sinai billiard in three dimensions proved beyond reasonable 
doubt that the methods developed for two dimensional chaotic systems can be 
extended to higher dimensions. Of utmost importance and interest is the study 
of classical chaos and its quantum implications in many body systems. This is 
probably the direction to which the research in "quantum chaos" will be advanc- 
ing. 
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A Efficient quantization of billiards: BIM vs. 
full diagonalization 

In this Appendix we wish to compare two possible quantization schemes for bil- 
liards: Direct Diagonalization (DD) of the Hamiltonian matrix vs. the Boundary 
Integral Method (BIM) (see e.g. ]T^, [88]). The diagonalization is a generic method 
to solve the time independent Schrodinger equation, while the BIM is specialized 
for billiards. To compare the two methods, we estimate the complexity of com- 
puting all of the eigenvalues up to a given wavenumber k. 

To find the matrix elements of the Hamiltonian we treat the billiard bound- 
aries as very high potential walls. The linear dimension M(k) of the Hamiltonian 
matrix that is needed for finding eigenvalues around k is: 



M m (k) = 0[ - =0((kS) d ) (185) 




where S is the typical linear dimension of the billiard, A = 2n jk is the wavelength 
and d is the dimensionality of the billiard. The above estimate is obtained by 
enclosing the billiard in a hypercube with edge S and counting the modes up to 
wavenumber k. The numerical effort to find eigenvalues of a matrix is of order 
of its linear dimension to the power 3. Thus, the numerical effort to find all the 
eigenvalues of the billiard up to k using DD is estimated as: 

C m (k) = 0{M* o {k)) = 0{(kS) M ). (186) 

The expected number of eigenvalues up to k is given to a good approximation by 
Weyl's law, which for billiards reads: 

N{k) = 0({kS) d ). (187) 

Thus, the numerical effort to calculate the first (lowest) N eigenvalues of a billiard 
in d dimension in the direct Hamiltonian diagonalization is 

C m (N) = 0(N 3 ) (188) 

which is independent of the dimension. 

As for the BIM, one traces the fc-axis and searches for eigenvalues rather than 
obtaining them by one diagonalization. This is done by discretizing a kernel 
function on the boundary of the billiard and looking for zeroes of the resulting 
determinant. The linear dimension of the BIM matrix is 

M Bm = ( (f ) j = ° (ikSy- 1 ) = O (N 1 - 1 /*) . (189) 
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This estimate is obtained from discretizing the boundary of the billiards which is 
of dimension d — 1 by hypercubes of edge A. The numerical effort of calculating 
the determinant once is: 

c B1M (k) = 0(Ml IM {k)) = 0((kS) 3 ^) . (190) 



(In practice, one often uses the SVD algorithm [ 55f , which is much more stable 
than a direct computation of the determinant and has the same complexity.) 
Using the relation ( |187f ) we find that the numerical effort to find an eigenvalue 
near the iVth one is estimated by 

c BIM (iV) = 0(N^ d ). (191) 

To get the above result we assumed that a fixed number of iterations (evaluations 
of the determinant) is needed to detect each eigenvalue, which is justified at least 
for the case where level repulsion is expected. Thus, the complexity to calculate 
all the eigenvalues up to the iVth is 

C Bm (N) = 0(N)c Bm (N) = 0{N A ~' i ' d ). (192) 

In particular: 

C B1M (N) = 0(N 5 / 2 ), forrf = 2 
C B1M (N) = 0(N 3 ), for d = 3. 

We conclude that the BIM is more efficient than DD for 2 dimensions, and 
for 3 dimensions they are of the same level of complexity. In practice, however, it 
seems that the BIM is better also in 3 dimensions, since the DD matrices can be 
prohibitly large, and manipulating them (if possible) can be very expensive due 
to memory limitations (paging). Also one has to take into account, that due to 
evanescent modes, the numerical proportionality factor in ( |189| ) is actually close 
to 1, while for ( |185| ) the factor can be large if high accuracy is desired. This is 
due to the fact that the off-diagonal matrix elements of the Hamiltonian decay 
only like a power-law due to the sharp potential and hence very large matrices 
are needed in order to obtain accurate eigenvalues. 
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B Symmetry reduction of the numerical effort 
in the quantization of billiards 

Consider a d-dimensional billiard which is invariant under a group Q of geomet- 
rical symmetry operations. We want to compare the numerical effort that is 
needed to compute the lowest N eigenvalues of the fully symmetric billiard with 
that of computing the lowest N eigenvalues of the desymmetrized billiard. In 
"desymmetrized" we mean the following: if Q is the full billiard domain, then 
the desymmetrized billiard u is such that Ugee 9 U = ^- ^ one uses ^ ne direct 
diagonalization (DD) of the Hamiltonian matrix, then there is no advantage to 
desymmetrization, because the prefactor in ( |188| ) should not depend on the shape 
of the billiard if its aspect ratio is close to 1. Therefore, the numerical effort of 
computing the lowest N levels of either the fully symmetric or the desymmetrized 
billiard is more or less the same using DD. On the other hand, as we show in 
the sequel, desymmetrization is very advantageous within the framework of the 
BIM,. 

We first note, that considering a particular irreducible representation 7 of Q is 
equivalent to desymmetrization of the billiard together with imposing boundary 
conditions that are prescribed by 7. The dimension of 7 is denoted as c? 7 and the 
order of Q is denoted as Ng. Given a complete basis of functions in which the 
functions are classified according to the irreps of Q, then the fraction of the basis 
functions that belong to the irrep 7 is d^/Ng = F 1 . This is also the fraction of 
eigenvalues that belong to 7 out of the total number of levels, when we consider 
a large number of levels. Using the notations of appendix [A], we thus have: 

Mfi&(A;) = F 7 M BIM (k) 

N^\k) = F y N(k) (193) 

where the quantities with superscript 7 correspond to the desymmetrized billiard, 
and the others to the fully symmetric one. Using (|189 ) and repeating arguments 
from appendix |A] results in: 

C£L{N) = FjC B1M (N) . (194) 

In the equation above we replaced — > N. Thus, the decrease in the density 

of states is more than compensated by the reduction in the size of the secular 

3 

matrix and the overall numerical effort is diminished by a factor of _F 7 d . For 
example, in the case of the 3D Sinai billiard and for a one-dimensional irrep, the 
saving factor is: 

1 /1 2 \^ 1 

F * = («) (195 > 

which is a very significant one. 
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C Resummation of Dim using the Ewald sum- 
mation technique 

In general, the Ewald summation technique is used to calculate (conditionally 
convergent) summations over lattices {p }: 



(196) 



One splits the sum S into two sums Si, 5*2 which depend on a parameter rj: 

s = Sl + s 2 = J2 h(p ; v) + E ^) • ( 19? ) 



This splitting is usually performed by representing f(x) as an integral, and split- 
ting the integral at 77. The idea is to resum 5*1 on the reciprocal lattice {g } using 
the Poisson summation formula: 



Si = ^2 J dd P ex p( 2ni P9)fi(p;v) = ^2h{g;v) , 



(198) 



and to choose 77 such that both Si and S 2 will rapidly converge. 

We need to apply the Ewald summation technique to D LM (k), given explicitly 
in equation @, which constitute the main computational load. This is because 
they include summations over the Z 3 lattice which need to be computed afresh 
for each new value of k. It is possible to apply the Ewald technique directly to 
each Dim, but it is much simpler to take an indirect route: We shall Ewald resum 
the free Green function on the 3-torus, and then read off the Z)la/'s as expansion 
coefficients. 

We start with the free outgoing Green function on the three-dimensional 
torus: 

yT 1 ^ exp(ik\q-p\) 



G » (?) = -£ E 



19- Pi 



where we took the side of the torus to be 1 for simplicity and defined q 



(199) 



To split the sum we use an integral representation of the summands |55|, p0| : 



exp(ik\q — p\ 

\q-p\ 



exp 



0(C) 



{q-pfe+ 



4£ 2 



(200) 



where the integration contour C is shown in figure 57. It is assumed that k has 
an infinitesimal positive imaginary part, which is taken to at the end of the 



calculation. We now deform the contour into C (see figure 17), such that it runs 
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along the real axis for £ > y/rj/2, and split the integral at this point as follows: 



cm 

G T 2 {q) 



Gl{q) + G T 2 {q) 



—T 



2ttv^F 



-y 



0(C) 
oo 



exp 



exp 



<q-p) 2 e+ 



4 £2 



d£, 



d£- 



(201) 
(202) 

(203) 
(204) 

The summation in is rapidly convergent, due to the fact that we integrate 
over the tails of a rapidly decaying function in £ (faster than a Gaussian), and 
we start further on the tail when p grows. In order to make G\ also rapidly 
convergent, we need to Poisson resum it. We use the identity: 

Er r- -a 2X2i n V^ ST^ r(27rr^ ! 
exp [-(q-p) £ J = —3- 2^ ex P 

p 9 



4£2 



i{2irg)q 



(205) 



which is obtained by explicitly performing the integrals of the Poisson summation. 
Thus, 



Gl{q) 



- J^exp(27ri#?) 



0(C) 



d£ 

^exp 



k 2 - (2ng)< 
4 £2 



exp(2iiigq) exp 



fc 2 -(27rg) 2 



A; 2 - (2vr^ 



(206) 



The second line was obtained from the first one by performing the integrals 
explicitly. The expression obtained for G\ is also rapidly convergent, and is 
suitable for computations. We thus succeeded in rewriting Gq as two rapidly 
converging sums ( |206| ), ( |203| ). We note that the results ( [206] ), ( |203| ) are valid for 
general lattices, the cubic lattice being a special case p0| . 

The heart of the above resummation of Gq was the integral representation 
( |200|) which is non-trivial. In appendix |D| we present an alternative derivation of 
the above results using more intuitive, physical arguments. 

It remains to extract the D^'s from the resummed Gq. The basic relation 



is 



k n (kq) 
D L M{k) + -7= . 71 N d L0 



4vr jo(kq) 



G T Q {q) = Y,JL{kq)Y LM {n 

LM L 

Using expansion theorems [5B| applied to ( [203| ) and 

G o(<?) = E y ^(<V) \ ^^ L e k2h Yt M {^)j L {2ixgq) ^ - ^ 



(207) 



one can rewrite as: 

-(2^) 2 A? 
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E -7=^(0?) / ^ d£ j L (-2 W £ 2 ) exp 

^ A/7T Jv^ 



W)£ 2 +^ 



fc 2 
4^2 



(20£ 



Comparing equations (p07|) and (|208|) , and using the orthogonality of the spherical 
harmonics Yla/(^V), one obtains: 



D LM {k) 



3L(kq) 



k 



n (kq)5 



LO 



(209) 



This is the Ewald-resumed expression of DiM(k). It has the interesting feature, 
that even though each of the terms explicitly depends on q, the total expression 
is independent of q. The same applies also to 77. This freedom can be used 
to simplify the expression ( |209| ), since for q — > the spherical Bessel functions 



simplify to powers p0|: 



him) 



(aqY 



(210) 



(2L + 1)!! ' 

which are computationally less demanding. Taking the limit is straightforward 
for L 7^ 0, while for D 00 there is a complication due to the singularity of n (kq). 
As shown in appendix [E] this singularity is exactly cancelled by the p = term, 
resulting in a finite expression also for Dqq. The final result is: 



D 



LM 



7-)(l) , T}W , T}W c 



(2) 



'LM 



LM 



'00 °i0 



(i) 

LM 



(2) 
LAf 



Am L k- L e k2 /^^9) L yt M {^ 
2 L+i k -L 



a 



(3) 
00 



, — 00 



(k 2 /r]) n 



^1 



.)— 

g, k 2 - (2ng) 2 
d£ i 2L exp 



Ve 2 + ^ 



4? 



2tt ^ n!(2n - 1) 

n=0 v 



(211) 
(212) 

(213) 
(214) 



with the convention g L \ g=0 L=0 = 1. This completes the task of Ewald-resuming 
the building blocks -Dlm(^) into rapidly convergent series. 
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D "Physical" Ewald summation of G^{q) 



In this appendix we present a derivation of the results ( |206| ), ( |203| ) by a method 
that is different than the one used in appendix 0. The present method is phys- 
ically appealing and does not require the use of complicated integral represen- 
tations. It is inspired by Appendix B of Kittel's book |97j which deals with 
the problem of calculating Madelung constants (electrostatic potentials) of ion 
crystals. 

In the sequel we use q = r — r ' and adopt the following notational convention: 
For any quantity X(q) we add a superscript T to denote its lattice sum: 

X T {q)^X{q-p). ( 215 ) 
p 

We start from the Helmholtz equation for Go: 

(Vl+k 2 )G (q)=5(q). (216) 
Due to linearity, the function Gq satisfies: 

(Vl+k 2 )G T (q)=6 T (q). (217) 



The RHS of ([217D can be interpreted as a "charge distribution" which is composed 



of point charges on a lattice. Each such point charge S(q—p) induces a "potential" 
Go(q-p) = — exp(ik\q— p\) / (4n\q— p\) which is long-ranged due to the sharpness 
of the charge. (This is in analogy to the electrostatic case.) Hence, the lattice 
sum of potentials Gq is conditionally convergent. To overcome this difficulty we 



introduce an arbitrary charge distribution X(q) and rewrite (|217|) as: 



Gl{q) = Gj(q)+G^(q), (218) 
(V 2 + k 2 )Gl (q) = X T (q), (219) 
(V 2 + k 2 )G T 2 {q) = 5 T (q)-X T (q). (220) 

We want X(q — p) to effectively screen the 5(q — p) charges, making G 2 short- 
ranged. This will result in rapid convergence of G\. (Note, that the equations 
( |218|) -( p20| ) hold also for the quantities without the T superscript due to linearity.) 
On the other hand, X(q) must be smooth enough, such that will rapidly 
converge when Poisson resumed. It is hence plausible to choose a (spherically 
symmetric) Gaussian charge distribution for X(q): 

X(q) = Aexp(-aq 2 ) , (221) 

where A and a are yet arbitrary parameters. 
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We calculate first G 2 (q) by rewriting the inducing charge as an integral over 
S charges, and using the fact that each 5 charge contributes G to the potential: 



(V 2 + k 2 )G 2 {q) = 8(q) - \(q) = 8(q) - / d 3 Q\(Q )8(q - Q ) . 



(222) 



Hence, 



G 2 (q) = G (q)- d 3 QX(Q)G (q-Q 



G (q) 



, 7l\- i k 2 

I - A ( -J 



+ 



A 
2aq J 



dt e~ a{t+q)2 cos(fct). (223) 



The first term is long-ranged due to Go, and the second term is short-ranged due 
to the integral that is rapidly decreasing as a function of q. To make G 2 short 
ranged, we thus have to set the coefficient of Go to 0, which is satisfied if we 
choose 



A = A(*,a)=(2)'axp(£). 

Therefore, we get for G\ a rapidly convergent sum: 

k 2 



(224) 



G T 2 {q) 



E 



27T0r ^ \q-p \ Jo 



dt exp [—a (t + \q-p\) ] cos(Jfci). (225) 



This can be re-expressed in a more compact form using complement error func- 
tions with complex arguments: 



Gl{q) 



where 



■-E 



1 



2tx * — ' \q — p\ 

p 



-Re 



exp(— ik\q — p\) ■ erfc ( y/a\q — p \ — 



ik \ 



1 f°° _ 2 
erfc(z) = —= / e u du . 

V 71 " Jz 



(226) 
(227) 



To calculate G\ we can directly Poisson resum (|225|) . Alternatively, we can 
use again the Helmholtz equation for G\ ( p!9|) to simplify the calculations. We 
expand G\ in the reciprocal lattice: 

G\{q) = J2j d3 P exp(2vr^)G 1 (g-p) 
9 

= ^2exp(2mqg) J d 3 p exp(-2vripp;)Gi(p ) 
9 

= J2 e M^Q9)Gi S (228) 
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where the second line was obtained from the first one by shifting the origin of 
the integration. Similarly for X T (q): 



(229) 



Inserting ( |228j , |229| ) into Q219| ) and using the orthogonality of the Fourier com- 
ponents, we get the simple relation between Gig and A^: 



G< 



ng k 2 -(2irg) 2 ' 
When inserted back into ( 228|) we finally get for G\ : 



(230) 



cm = e 



exp(27cigq) exp 



fc 2 -(27rg) 2 
4a 



k 2 - (27rg) 2 



(231) 



This expression is identical with ( |206| ) if we set 4a = r\. It can be shown 
that also the expressions for G\, ( p03| ) and ( |226| ) are identical. However, equa- 
^2^) is more convenient if one needs to compute G^(q), since it involves 



tion 

well-tabulated computer-library functions |58[ and saves the burden of numerical 
integrations. On the other hand, the expression fl203p is more convenient as a 
starting point for calculating D LM {k). 

To summarize, we re-derived the Ewald-resummed form of G^(q) using the 
underlying Helmholtz equation. We used a physically intuitive argument of 
screening potentials, that was shown to be equivalent to the more abstract inte- 
gral representation of G (q), equation (|20( 
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E Calculating D t 



(3) 
00 



We need to calculate (refer to equation ( |209|) and its subsequent paragraph): 



= lim 



3o(kq) 



1 cos(A;g) 1 

7T 



47r q 



(232) 



where we used the explicit expression no(x) = —cos(x)/x. Taking the limit of 
the denominator is trivial, since ja(kq) — > 1. For q — > we can write, 



1 cos(kq) 



47r q 



Arcq 



+ 0{q) 



(233) 



which contains \jq singularity. As for the term with the integral, we expand 
exp(/c 2 /4£ 2 ) in a Taylor series, and transforming to the variable t = q£ one gets: 



- / d^exp - g 2 e 2 + — = — J2^rr \ 



dt t"^e 



2n -t 2 



(234) 



For n = 0: 



00 riV^/ 2 



dt e~* = 
For n > we integrate by parts: 



dt e" 



v% + 0(g 2 ). (235) 



dt t~ zn e 



In 



-2n+l 



-2n+3\ 



(236) 



Collecting everything together back to (|232|) , the 1/q singularities cancel, and we 
remain with the finite expression: 



a 



(3) 

oo 



2tt ^ n!(2n - 1) 

n=0 v y 



(237) 
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F The "cubic harmonics" Yj^j K 

F.l Calculation of the transformation coefficients ai) KM 

We want to calculate the linear combinations of spherical harmonics that trans- 
form according to the irreducible representations of the cubic group Oh- This 



problem was addressed by von der Lage and Bethe |57j which coined the term 
"cubic harmonics" for these combinations. They gave an intuitive scheme that 
was used to calculate the first few cubic harmonics, but their arguments are dif- 
ficult to extend for large Us. Moreover, their method is recursive, because one 
has to orthogonalize with respect to all lower lying combinations. This is cum- 
bersome to implement numerically and might result in instabilities for large L's. 
The only other work on the subject that we were are aware of specializes in 
the symmetric representation and gives only part of the combinations. It also ex- 
presses the results not in terms of spherical harmonics, but rather as polynomials 
that are difficult to translate to Ylm's. 

We describe in the following a simple and general numerical method to cal- 
culate the cubic harmonics in a non-recursive way. This is based on a general 
theorem, that states that a function transforms according to the irrep 7 iff 



it satisfies ||31|| : 



= fbt) (238) 
where P^) is the projection operator onto the subspace that belongs to 7: 

We denoted by / 7 the dimension of 7, Nq is the number of elements in the group 
G, and x id) are the characters. The realization of P^ as a matrix in an 
arbitrary basis results in general in an infinite matrix. However, in the case of 
the cubic harmonics, we know that Oh C 0(3), thus the operations of g e Oh do 
not mix different L's. Hence, working in the Ylm basis, we can write the cubic 
harmonics as the finite combinations: 

Y£>m = £ a%Y LM (n) (240) 

M=-L 

where J enumerates the irreps 7 in L. For simplicity we consider 1-dimensional 
irreps. Applying ( [238] ), ( 239 ) to (|240|) and using the Wigner matrices VW(g) to 



express the operations of g on Ylm pEfl , we get the following (2L + 1) x (2L + 1) 
linear system: 

" "2m- = (241) 



E 

M' 



MM' ~ "MM' 
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where: 

«4E^(#«(^ (242) 

g&G 

The above equations are best solved using SVD algorithm and the (orthonor- 
malized) eigenvectors that belong to the zero singular values are the required 
coefficients o^)* M ,- For multi-dimensional irreps one needs to classify the cubic 
harmonics also with respect to the row K inside the irrep. This can be done 
by simple modification of the above procedure, using the appropriate projectors 

The above general procedure can be simplified for specific irreps. In the 
following we shall concentrate on the completely antisymmetric irrep 7 = and 
further reduce the linear system ( [24 We first note, that the antisymmetric 
cubic harmonics must satisfy per definition: 

gY&\n) = X (a) (g)YS(ty = (-l)^° !S) Yff(Q) Vg e O h . (243) 

We then choose a few particular g's for which the operations on Ilm(^) ar e 
simple: 

f x (xyz) = (-xyz) : f x Y LM (6, 0) = Y LM (6, -0) = (-l) M Y L _ M (d, 0) 

f y (xyz) = (x- yz) : f y Y LM (6, 0) = Y LM (6, n - 0) = Y L _ M (6, 0) 

f z (xyz) = (xy - z) : f z Y LM (6, 0) = Y LM (n - 9, 0) = (-1) L+M Y LM (8, 0) ^ V 



p xy (xyz) ee (yxz) : p^Y LM (9, 0) = Y LM (9, | - 0) = (-i) M Y; 



L-M{V, > 



Applying ( |243| ) and ( |244| ) to (|240D results in the following "selection rules": 

«Sm = 0, L^2p+l,M^4q, p,qeJN 
aif-M = (245) 

which reduces the number of independent coefficients to be computed by a factor 
of 16. The form of the projector matrix p( ai ) can also be greatly reduced, if we 
observe that the group Oh can be written as the following direct multiplication: 

O h = G 3 ®G 16 (246) 
G 3 = {e,c, c 2 } e = identity, c{xyz) = (yzx) (247) 
Gie = {e, Pxy] ® {e, r x }<8 {e, f y } <g> {e, f z } (248) 

and consequently, the projector can be written as 

P<°> = AAe (249) 
p 3 = e + c + c 2 (250) 
Ae = {e-p xy ){e-r x )(e-f y )(e-f z ). (251) 
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The operator Pi 6 acts as the identity on the subspace defined by ( |245| ) and hence 
we need to consider only the operation of P3. Simple manipulations give the 
following set of equations: 

E [ 2d SW (£) - <W] a[% ql =0, g > , L — 2p+l. (252) 

<?'>0 

The matrices , are the "reduced" Wigner matrices, which are real |p6[ , thus 
the resulting coefficients are also real. The above is a square linear system, which 
is 8 times smaller than the general one ( 241|) . 

F.2 Counting the Y^'s 

The number of the irreps 7 of Oh that are contained in the irrep L of 0(3) is 
given by the formula [ |31|j : 

N,L = ^J2^(9)XL(g) (253) 

geO h 

where Xl{q) are the characters of the irrep L. An explicit calculation shows that 
the main contributions for large L's come from the identity and from the inversion 
operations, thus: 

N, L n[l±{-lf] M2 ^ + 1) . (254) 

where the ± corresponds to the parity of 7. Since for Z 7 -dim irrep we have Z 7 
basis functions, and there are 2L + 1 basis functions in the irrep L, the fraction 
of cubic harmonics that belong to 7 is: 

F y « ^ (255) 
in accordance with the general relation: 

E/? = 48. (256) 

7 

Consequently, the fraction of cubic harmonics that belong to the i^'th block of 
7 is 

F,k « i ■ (257) 
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G Evaluation of l(p p ) 

G.l Proof of equation ( |ID|) 

We need to prove the relation: 



Y, Y lm{^p) = Kp p ) *W*V) (258) 

9&O h p'eS(pp) 

where p p = k) is the unique vector in the set OhP which resides in the fun- 
damental domain % > j > k > 0, S(p p ) is the collection of all distinct vectors 
obtained by the operations gp p , g G O^, and l(p p ) is an integer. 

Proof. Let TC be the set of all g G Oh under which p p is invariant: 

9P P = Pp 9 e H. (259) 

The set TC is a subgroup since: 

1. The identity e G TC. 

2. The set TC is closed under multiplication, since if g\, §2 G TC then gi(g2P P ) = 
9iPp = Pp. 

3. The set TC is closed under inversion: g~ l p p = g~ 1 (gp" p ) = p p - 

The order of (number of terms in) TC is denoted as iV%. We assume that 7^ is the 
maximal invariance subgroup, and construct the right cosets gTC = {gh±, . . .}. Ac- 
cording to |3T| there are N c = 48/iV^ mutually exclusive such cosets Ci, . . . , Cn c 
(The number 48 is the order of Oh). Their union is Oh- For each coset Cj we can 
define 

Pi = dp p (260) 

which is meaningful due to the invariance of p p under TC. 
We want to prove the following 
Lemma, pi 7^ pj iff i ^ j. 

Proof. Assume the opposite, then in particular 

9iP v = 9jPp 

O (9j l 9i)p P = P P 

& gj 1 ^ = heTC 

9i = 9jh 
4=^ C[ = Cj 
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in contradiction to the assumption. The last line was obtained using the rear- 
rangement theorem applied to the group H. QED. 
We now set 



S(p p ) = (261) 
l(p p ) = N n = integer (262) 

and since Y,g & o h = Efii E^q we proved (|5|). QED. 



G.2 Calculating /(/%,) 

We give an explicit expression of l(p p ). Consider p p = (i,j, k) such that i > j > 
k > with no loss of generality. Then 

KPv) = lp(Pp)ls(Pp) , (263) 

( 1, i^j^k^i 
1 p(Pp) = \ 2, i = j ^ k or % ^ j = k or % = k ^ j (264) 
[6, i= j = k 

l s (p p ) = 2 (# zeroindices ). (265) 

We prove this formula in the following. First we observe, that Oh can be decom- 
posed as 

O h = V 3 ®S 3 (266) 
P 3 = group of permutation of 3 numbers (267) 
<S>3 = {± ± ±} = 3 sign changes. (268) 

Let Hp,Hs be the subgroups of V3,Ss, respectively, under which p p is invariant. 
Lemma. H = Hp ® Hs- 

Proof. Let g = ps, where g G Oh, p G V3 and s G S3. This representation 
of g is always possible according to Q266|) . If s $(Hs then sp p 7^ p" p , thus neces- 



sarily there is at least one sign change in sp~p with respect to p p . Consequently, 
gp p 7^ p~p, because permutations only change the order of indices and cannot re- 
store the different sign(s). We conclude that g ^H. Thus, g <EH => s E Hs- For 
every g G H we must have therefore gp p = psp p = pp p = p p which proves that 
also p G Hp. QED. 

We conclude that Nn = order (Hp) ■ order (Hs)- This is manifest in equations 
(|263| - g65| ). 
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H Number theoretical degeneracy of the cubic 
lattice 



H.l First moment 

The following arguments are due to J. Keating |9l| . We first need to estimate 
the fraction of integers that can be expressed as a sum of 3 squares. The key 
theorem is due to Gauss and Legendre and states that: 

q = f + f + k 2 , i,j,keJN g^4 m (8/ + 7), m,l e¥S. (269) 

From this we can estimate that the fraction of integers which can not be expressed 
as a sum of 3 squares of integers is: 

i( 1 + i + i. + ...) = i. (270 ) 
8 V 4 4 2 J 6 V ' 

In the above we used the fact (which is easily proven) that if q = 4 m (8Z + 7) then 
m, I are uniquely determined. Therefore, asymptotically only 5/6 of the integers 
are expressible as a sum of three squares. 

Our object of interest is the degeneracy factor g p (p) defined as: 

g p { P ) = ^e^\K = p) . (271) 

The number of Z 3 -lattice points whose distance from the origin is between p and 
p + Ap is estimated by considering the volume of the corresponding spherical 
shell: 

N p « 4vrp 2 Ap . (272) 
Since p 2 is an integer, the number of integers in the same interval is: 

n p « 2pAp . (273) 

Taking into account that only 5/6 of the integers are accessible, we obtain: 

(9 P (P)) = = ■ (274) 



(5/6)™ 



p 



H.2 Second moment 

Here we use a result due to Bleher and Dyson |[LO0|l , brought to our attention by 
Z. Rudnick: 

N lQn 2 C(2) 

S2g 2 p (yJk) = cN 2 + error , c = ^^77^ ~ 30.8706 . (275) 
k=i 7 
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Differentiating by N and considering only integers for which g p ^ one obtains: 



Therefore, 



12 

<<?») ~ ^cp 2 » 74.0894p 2 . (276) 



^4^-~/9p, /? = - ^ 9.8264. (277) 
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I Weyl's law 



A very important tool in the investigation of eigenvalues is the smooth counting 
function, known as Weyl's law. For billiards it was thoroughly discussed e.g. by 
Balian and Bloch p2| and by Baltes and Hilf |9l | . We construct in the following 
the expression for the 3D Sinai billiard. In general, it takes on the form: 



N(k) = N 3 k 3 + N 2 k 2 + + N 



(278) 



where we included terms up to and including the constant term. In fact, for 
the nearest-neighbour and two-point spectral statistics the constant term Nq 
is unimportant, since it shifts the unfolded spectrum x n = N(k n ) uniformly. 
Nevertheless, for completeness we shall calculate this term. We enumerate the 
contributions in the case of Dirichlet boundary conditions one by one and then 
write down the full expression. Figure ^| should be consulted for the geometry of 
the billiard. 



N3: There is only one contribution due to the volume of the billiard: 

volume 1 



Q-r 2 



288vr 2 



S A 



-nR 3 
3 



(279) 



N2: The contribution is due to the surface area of the planes + sphere: 

surface 1 



No 



16n 



384tt 



6(1 + V2)S 2 - 7nR 2 



(280) 



Ni: Here we have contributions due to the curvature of the sphere and due to 
2-surface edges: 



curvature: 



NI 



curvature 



127T 2 



ds 



surface 



1 1 

+ 



Ri(s) B*(s) 



R 



(281) 



where R±, R2 are the principal local radii of curvature. 

edges: We have 6 plane-plane edges and 3 plane-sphere edges. Their con- 
tributions are given by: 



N{ 



edges 



24vr 



edges 



7T 



a, 



71 



(282) 



' S (27 + 9V2 + 8^) L 



1447T 



24vr V 8 



9tt 



95 
12 



where Lj are the lengths of the edges, and ay are the corresponding 
angles. 
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There are three terms here due to square of the curvatures, 3-surface corners 
and curvature of the edges: 



curvature 2 : 



N, 



curvature 2 







5127T 



ds 



surface 



Ri(s) 



(283) 



3-surface corners: In the 3D Sinai billiard we have 6 corners due to inter- 
section of 3 surfaces; 3 of them are due to intersection of 3 symmetry 
planes and the other 3 are due to intersection of 2 symmetry planes 
and the sphere. The corners are divided into 4 types as follows: 

1 x a = (45°, 54.74°, 36.26°) 

3 x (3 = (45°, 90°, 90°) 

1 x 7 ee (60°, 90°, 90°) 

1 x 5 ee (90°, 90°, 90°) . 



As for the corners /3, 7, 8 which are of the type 
known expression for their contribution 



C0 



1 /tt 
96 V0 



7T 



, 90°, 90°) there is a 



;284) 



Therefore: 



(285) 



5 1 1 

"128' C7 ~"36' C<5 ~~64- 
As for the corner a, we calculate its contribution from the R = 
integrable case ("the pyramid"). The constant term in the case of the 
pyramid is —5/16 [|53[] and originates only from 3-plane contributions 
(there are no curved surfaces or curved edges in the pyramid). The 
pyramid has 4 corners: 2 of type a and 2 of type f3. Using c@ above 
we can therefore eliminate c„: 



2 • c a + 2 



128 



5 

16 



15 



128 



(286) 



Hence, the overall contribution due corners in the 3D Sinai is: 

15 \ / 5 



3— surface 




1 



128 



+ 3 



128 



+ 



1 36; V 64 

5 

18 ' 



;287) 
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curvature of edges: We have 3 edges which are curved. They are 90° 
edges that are due to plane-sphere intersections. Baltes and Hilf [|9~1~| 



quote the constant term (— l/12) + (l/256)(if/-R) for the circular cylin- 
der, where H is the height and R is the radius of the cylinder. We 
conclude from this that the if-independent term —1/12 is due to the 
curvature of the 90° edges between the 2 bases and the tube. Assum- 
ing locality, it is then plausible to conjecture that the contribution due 
to the curvature of a 90° edge is: 



1 



48vr 



edge 



dl 

W) 



(2f 



where R(l) is the local curvature radius of the edge. When applied to 
our case (R(l) = —R), we get: 



Arcurv. edge 
iV 



1 

64 



(289) 



Putting everything together we get: 
N(k) -- 



— ^ [S 3 - -irR 3 ) k 3 
288tt 2 V 3 



+ 



384tt 

144vr 
151 
576 ' 



6(1 + v^S 2 - 771R 2 



(-- 


11 )] 




32vry_ 



(290) 
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J Calculation of the monodromy matrix 



The monodromy matrix measures the linear response to infinitesimal displace- 
ments of the initial conditions of a classical orbit. Its eigenvalues determine the 
stability of the orbit. Due to the symplectic form of the equations of motion, 
if A is an eigenvalue of the monodromy matrix then also A*, 1/A and 1/A* 0. 
Therefore, generically the eigenvalues come in groups of four: 

A = exp(±u ± iv) , u,veTR. (291) 

In d dimensions the monodromy has 2(d — 1) eigenvalues. Therefore, only for 
d > 3 the generic situation (|291|) can take place. In two dimensions there are 
only two eigenvalues and consequently one obtains the following three possible 
situations (which are special cases of ( [2 9 If ) with either u or v set to 0): 

Elliptic: Ai,2 = exp(±ii>), stable orbit. 

Parabolic: \i 2 = 1 or Ai^ = (— 1), neutrally stable orbit. 

Hyperbolic: Ai^ = exp(±w) or Ai 5 2 = — exp(±w), unstable orbit. 

The parabolic case with the "+" sign is denoted as "direct parabolic" and with 
"-" sign it is denoted as "inverse parabolic" . Similar terminology applies to the 
hyperbolic case. The generic case ( |291| ) is designated as "loxodromic stability" 

0- 



J.l The 3D Sinai torus case 

We wish to calculate explicitly the 4x4 monodromy matrices in the case of the 
3D Sinai torus. There are (at least) two possible ways to tackle this problem. 
One possibility is to describe the classical motion by a discrete (Hamiltonian, 
area-preserving) mapping between consecutive reflections from the spheres. The 
mapping is generated by the straight segment that connects the two reflection 
points, and the monodromy can be explicitly calculated from the second deriva- 
tives of the generating function. This straightforward calculation was performed 
for the 2D case (for general billiards) e.g. in [[72] and it becomes very cumbersome 
for three dimensions. Rather, we take the alternative view of describing the clas- 



sical motion as a continuous flow in time, as was done e.g. by Sieber [BU| for the 
case of the 2D hyperbola billiard. We separate the motion into the sections of free 
propagation between spheres and reflections off the spheres, and the monodromy 
matrix takes the general form: 

M = M^M^ ■ ■ ■ M^M^M^M^ , (292) 

where M^p*"' describes the free propagation from sphere i to sphere i + 1 and 
M* ef describes the reflection from the sphere i. To explicitly calculate the matrices 
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one has to choose a well-defined (and convenient) coordinate system, which is a 
non-trivial task in three dimensions. If we denote the direction along the orbit by 
"1", then we have two more directions, denoted henceforth "2" and "3". Hence 
there is a rotation freedom in choosing the directions 2 and 3. For convenience of 
calculation of M rei we choose the following local convention for coordinates: Near 
sphere % there exists the plane Vi which is uniquely defined (except for normal 
incidence) by the incoming segment, the outgoing segment and the normal to 
the sphere % at the reflection point. Direction 1 is obviously in Vi. We uniquely 
define direction 3 to be perpendicular to V% along the direction of the cross 
product of the outgoing direction with the normal. Direction 2 is then uniquely 
defined as e 2 = e 3 x e.\ such that a right-handed system is formed. Obviously 
e-2 is contained in Vi. The uniqueness of the local coordinate system guarantees 
that the neighbourhoods of the initial and the final points of the periodic orbits 
are correctly related to each other. To account for the local coordinate systems 
we need to apply a rotation between every two reflections that aligns the "old" 
system to the "new" one. Hence, 



M 



M n+1 ^ n M n t 1 ^ n M n f 

prop rot ret 



M 3^2 M 3<-2 M 2 M 2^1 M 2<-1 M 1 
prop rot ref prop rot ret 



(293) 



We should also fix the convention of the rows and columns of M in order to be 
able to write explicit expressions. It is chosen to be: 



/ Sq 2 \ 

$ <?3 
V <^3 J 



M 



final 



/ 5q 2 \ 

<5<?3 

V &P3 J initial 



(294) 



A detailed calculations gives the explicit expressions for M prop , M ref and M rot : 



lv± prop 



K* = 



( 


1 




1 


-i/p 












\ 












1 


£j+l<- 


-i/p 




\ 













1 




J 


/ 




-1 










°\ 








2p 


-1 
















R cos (3i 


















1 







V 










2p 


COS (3i 

R 


1 J 





cos ctj+i^j 

cosctj+i^ 

sin a i+1 ^i 

sin a i+ i<. 



- sin a i+1 ^i 


COSttj+i^i 







- sin a i+1 ^i 


COSftj+i^i J 



(295) 



(296) 



.(297) 



In the above p is the absolute value of the momentum which is a constant, L i+1 ^i 
is the length of the orbit's segment between spheres i and i + 1, fa is the reflection 
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angle with respect to the normal of the sphere i and is the angle that is 

needed to re-align the coordinate system from sphere i to i + 1. Even thought 
the entries of M are dimensional, the eigenvalues of M are dimensionless. Hence, 
the eigenvalues cannot depend on p, which is the only variable with dimensions 
of a momentum. (All other variables have either dimension of length or are 
dimensionless.) Therefore, one can set p — 1 for the sake of the calculations 
of the eigenvalues of M. The formulas above for the monodromy were verified 
numerically for a few gainst a direct integration of the equations of motion 

near a periodic orbit of the Sinai torus. We mention the work of Sieber [52| which 
extends the calculation of the monodromy matrix to an arbitrary billiard in three 
dimensions. 



J. 2 The 3D Sinai billiard case 



We next deal with the calculation of the monodromy matrix for the periodic orbits 
of the desymmetrized 3D Sinai Billiard. In principle, one can follow the same 
procedure as above, and calculate the monodromy as for the Sinai torus case, 
this time taking into account the presence of the symmetry planes. A reflection 
with a symmetry plane is described by: 



/1/fPlane 
JW ref 



/ 


-1 








o\ 







-l 
















1 





V 











1 / 
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which is simply M re f with R —>■ oo. This method, however, is computationally 
very cumbersome because of the need to fold the orbit into the desymmetrized 
Sinai billiard. Instead, we can use the monodromy matrix that is calculated for 
the unfolded periodic orbit, because the initial and final (phase space) neighbour- 
hoods are the same modulo g. A calculation shows, that in order to align the 
axes correctly, one needs to reverse direction 3 if g is not a pure rotation: 



w 



( 1 








o \ 





1 






























unfolded 

w 



(299) 



where a(g) is the parity of g: 



-1. 



g is a rotation 

g is an improper rotation (rotation + inversion) 



(300) 



The above formulas were verified numerically for a few cases by comparing the 
result ( |299| ) to a direct integration of the classical dynamics in the desymmetrized 
Sinai billiard. 
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